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I. INTRODUCTION 



I. 


INTRODUCTION 


The global atmospheric and oceanographic research activities 
of the Goddard Institute for Space Studies were transferred to the 
Modeling and Simulation Facility of the new Goddard Laboratory for 
Atmospheric Sciences (GLAS) in 1977. 

This Research Review presents a synopsis of studies conducted 
at the Facility during the past year. The Review is organized according 
to three research areas: global weather, ocean/air interactions, and 
climate. A list of recent publications by the staff scientists and 
visiting scientists is presented at the end of the Review. 
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II. GLO.BAL WEATHER 


$^C£Dm PAGE &AW P«T P2^. 



MODEL-DEPENDENCE STUDIES OF THE EFFECTS OF SOUNDING TEMPERATURE 
ASSIMILATION ON WEATHER FORECASTING 

{M. Ghil, M, Halem, and R. Atlas) 

In comparing various methods for the assimilation of 
remote sounding information into numerical weather prediction 
models, the problem of model dependence for the different results 
obtained is clearly important. Most assimilation experiments 
have not addressed this problem directly, every research and 
development group working essentially with one numerical model. 

GLAS has started to investigate two aspects of the model 
dependence question; (1) the effect of increasing horizontal 
resolution within a given model on the assimilation of sounding 
data, and (2) the effect of using two entirely different models 
with the same assimilation method and sounding data. In order to 
study the first aspect, GLAS has applied its four-dimensional 
(4-D) time-continuous statistical assimilation method (SAM, Ghil 
et ai . , 1979) to the assimilation of sounding data from the Data 
Systems Tests (DST) period DST-6 (February 1976) into the GLAS 
general circulation model with a coarse resolution (4°lat.x5°long. ) 
and with a fine resolution (2. 5°lat .x3°long. ) . To study the second 
aspect of model dependence, an asynoptic, time-continuous successive 
correction method (SCM, Ghil et al . , 1979) was applied to the global, 
ninelevel, primitive-equation model of the U. S. National Meteoro- 
logical Center (NMC, Stackpole, 1976). Results of SCM assimilations 
with the latter model are compared with those of previous SCM 
experiments with DST-6 data using the GLAS coarse-mesh model. 

Preliminary results^ of comparisons between SAM 
assimilations with the coarse-mesh and the fine-mesh GLAS model 
show that increased model resolution does not affect the impact 
of satellite data on weather forecasts in any negative way; the 
effects of increased resolution and of satellite data assimilation 
are cumulative: they both tend to improve forecasting accuracy by 

a similar fractional amount. This fact is reflected in the usual 
numerical measures of forecast accuracy, viz., RMS errors and Si 
skill scores. It is also reflected in synoptic-case studies of 
the large-scale forecast fields'. 

Assimilation cycles were run using data from January 29 
until February 21, 1976, with both the coarse and the fine GLAS 
models. For each model resolution, one cycle was run including 
all data from the NOAA-4 and Nimbus-6 satellites, and another 
one was run with the same data, but excluding the temperature 
sounding information: they are referred to, respectively, as 

the 2SAT and NOSAT cycle. For each one of these four experiments, 
72-hour forecasts were performed from the analyses produced by 
the corresponding assimilation cycle at 48-hour intervals, 
starting on February 1. The forecasts were verified at 48 hours 


3 



and at 72 hours against NMC operational analyses. The quantities 
verified and regions of verification are the same as in 
Ghil et al . (1979) . 

The results are summarized in Table 1, according to 
the methods outlined in Ghil et al (1979). It is seen that the 
effect of model improvement, specifically increased resolution 
on the one hand, and utilization of satellite data on the other, 
act both in the direction of improving the forecast accuracy. In 
fact, the two effects appear to be nearly additive, as indicated 
by a comparison of the results in the last three rows of the 
table. Note also that the 2SAT experiment 2F, with the fine 
model, consistently gave impacts that are significantly higher 
than those of the corresponding experiment with the coarse 
model 2C, when measured against the common standard of the 
coarse NOSAT experiment OC (viz. the second and the last row 
of the table) . 


Table 1. Summary of Results for Model-Dependence Experiments. 


Experiment 

NO. 

Comparison 

Code 

Impact 

Si 

(%) . 
RMS 

Statistical Signif. 
(average/standard error) 

Sj_ RMS 

OC 

OC 

0 

0 



2C 

2CC 

4.46 

12.41 

1.85 

2.. 5 9 

OF 

OF 

0 

0 

_ 



OFC 

5.86 

7.30 

1.62 • 

1.35 

2F 

2FF 

4.49 

7.98 

1.40 

1.71 


2FC 

10.01 

14.47 

2.67 

2.47 


The symbol 0 stands for a NOSAT experiment, 2 stands for 
a 2SAT experiment; the symbol C stands for an experiment with the 
coarse model, F for one with the fine model. The results of the 
2SAT experiment with the fine model 2F were compared with those 
of the coarse NOSAT experiment OC and with those of the fine 
NOSAT experiment OF; these comparisons are denoted by 2FC and by 
2FF, respectively. The NOSAT experiment OP appears both as a 
control for 2F, viz,, as OP, and in a comparison versus the coarse 
NOSAT OC, viz., as OFC. 
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Previous results with the'GLAS coarse-mesh model 
(Ghil et al » / 1979) indicated that by using a time-continuous 
4-D method with statistical features (SAM or SCM) the assimi- 
lation of remote sounding temperature data can lead to numerical 
forecasts that show a modest, but statistically significant, 
improvement over forecasts started from analyzed fields not 
containing the remote sounding information. To determine 
whether the translation of GLAS methods into operational 
practice is feasible, it was decided to apply them to the 
model used by NMC to produce their operational analyses 
(Stackpole, 1976), which shall be called in the sequel the 
9LG (nine-level global) model. 

As a first step, SCM was adapted to and reprogrammed 
for the 9LG model. Numerous problems arose in the adaptation, 
due to differences between the models, and some of these problems 
await better solutions than the ones given to them at the present. 
NOSAT and 2SAT assimilation cycles were started and carried out 
from February 1 until February 7, 1976, using the same data sets 
as for the corresponding GLAS model cycles. Forecasts were 
carried out from the analyzed fields generated on February 3, 

5, and 7. These experiments are being continued and will be 
discussed more completely in a future publication. 

At this point, a comparison of results with those of 
SCM experiments performed on the GLAS model (Ghil et al . , 1979) 
shows that the differences between analyzed fields from the 
NOSAT and 2SAT cycles do not seem to' depend on the model: the 

differences produced by the utilization method of satellite 
data are of similar magnitude for the 9LG model as they were 
for the GLAS 4°x5® model. These differences range up to 100 ra 
in the 500 mb geopotential height field, especially over the 
North Pacific, where conventional data coverage is relatively 
poor. These initial state differences lead also to correspond- 
ingly large differences in 48-hour and 72-hour forecasts. GLAS 
intends to examine further the question of whether or not these 
differences are due solely to the assimilation of satellite data. 

The GLAS experiments on model dependence of sounding 
temperature assimilation suggest the following tentative 
conclusions : 

1. That model improvement, as exemplified by increased 
resolution, can act in the same direction as judicious 4-D assimi- 
lation of remote sounding information to improve 2-3 day numerical 
weather forecasts. The effect of model improvements and the 
effect of an observational data base augmented by the proper 
inclusion of sounding temperatures, when combined, can lead to 
numerical prognoses which achieve improvements of meteorological 
significance more easily than either effect separately, 

2. That the time-continuous 4-D methods developed at 
GLAS have similar beneficial effects when used in the assimilation 
of remote sounding information into NWP models with very different 


5 



numerical and physical characteristics. This suggests the 
possibility of the eventual implementation of these time- 
continuous methods in an operational environment, in which 
they might contribute to noticeable improvement in routine 
weather forecasting. 

References 


Ghil, M., R. Billing, and H. Carus , 1977: Proc. 5th Conf. 

Probability Stat. Atmos. Sci ., Amer. Meteor. Soc., Boston, 

MA, pp. 320-324. 

, M. Halem, and R. Atlas, 1979: Time-continuous assimilation 

of remote- sounding data and its effect on weather forecasting. 
Mon. Wea. Rev ., in press. 

Halem, M., M. Ghil, R. Atlas, J. Susskind, and W. J. Quirk, 1978: 

The GISS Sounding Temperature Impact Test . NASA Tech . Memo . 
78063, 421 + viii pp. 

, M. Ghil, and R. Atlas, 1978: Some experiments on the 

effect of remote sounding data upon weather forecasting. 

Remote Sensing of the Atmosphere: Inversion Methods and 

Applications , A. L. Fymat and V. E. Zuev., eds., Elsevier, 
Amsterdam. ■ 

Stackpole, J. D. , 1976: Proc. 6th Conf. Weather Forecasting Analysis , 

Amer. Meteor. Soc., Boston, MA, pp. 112-116. 
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SYNOPTIC EVALUATION OF THE IMPACT OF SOUNDING DATA INSERTIONS 
AND INCREASED HORIZONTAL RESOLUTION ON GLAS MODEL FORECASTS 

{R. Atlas, R. Rosenberg, and M. Eaton) 


The inpact of satellite-derived temperature soundings 
on short- tern (0-72 hour) weather forecasting has been expressed 
in terns of statistical measures of .forecast accuracy, regional 
differences in prognostic flow patterns, and automated local 
forecasts of precipitation occurrence (Halera, et al . , 1978). As 
an additional measure of the impact on forecast utility, a synoptic 
evaluation was performed by directly comparing the prognostic charts 
generated from 2SAT and NOSAT initial conditions (with and without 
satellite data insertions respectively) with each other and with 
the corresponding NMC analyses for verification. Comparisons were 
made every 12 hours for a variety of prognostic fields for eleven 
different forecasts from February 1976. The results of our com- 
parisons for the sea-level pressure and 500 mb height fields 
generated from the GLAS 4° latitude by 5® longitude model and 
verified over North America are presented in Table 1. 

This table gives the consensus of experienced forecasters 
in their subjective assessment of the prognostic charts. During 
this phase of the evaluation, the forecasters had no knowledge of 
which prognostic charts were based on 2SAT initial conditions. 
Following Atkins and Jones (1975), each forecast has been classified 
according to the following scale: 

A - 2SAT significantly better than NOSAT, 

B - 2SAT better than NOSAT, 

C - 2SAT and NOSAT of equal quality, 

D - 2SAT worse than NOSAT, 

E - 2SAT much worse than NOSAT. 

This classification system is designed to represent the 
relative utility of the prognostic charts generated from 2SAT and 
NOSAT initial conditions. The term significant in categories A 
and E implies that major prognostic differences have occurred. 
Categories B and D imply that there has been only a minor alteration 
of the prognostic pattern or that the net impact for the entire 
forecast region is small. 

Prom the table it can be seen that at sea level, none of 
the prognostic charts have been classified as A or E and that an 
almost equal number of B and D marks have been assigned. Although 
no significant impact has occurred at this atmospheric level, a 
tendency for more beneficial impacts at 48, 60, and 72 hours is 
evident . 


At 500 mb, two 72-hour prognostic charts have been 
classified as A and there were no E marks. In addition more than 
three times as many B than D marks have been assigned. Impacts 
at this level are most consistently beneficial at 60 and 72 hours. 
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Table 1. Results of a subjective assessment of 2SAT and 
NOSAT forecasts using the GLAS 4‘’x5° model. 


la. Sea Level Pressure 


CATEGORY 

12h 

24h 

36h 

48h 

60h 

72h 

TOTAL* 

A 

0 

0 

0 

0 

■0 

0 

0 

B 

0 

1 

1 

4 

3 

4 

13 

C 

8 

8 

7 

5 

6 

5 

39 

D 

3 

2 

3 

2 

2 

2 

14 

E 

0 

0 

0 

0 

0 

0 

0 


lb. 

500 mb 

Geopotential Height 


1 

CATEGORY 

12h 

24h 

36h 

48h 

60h 

72h 

TOTAL* 

A 

0 

0 

0 

0 

0 

2 

2 

B 

2 

2 

2 

3 

4 

4 

17 

C 

8 

9 

9 

5 

7 

4 

42 

D 

1 

0 

0 

3 

0 

1 

5 

E 

0 

0 

0 

0 

0 

0 

0 


* Total number of prognostic charts in each category. 
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A synoptic evaluation in the sane manner as that 
described above has also been conducted in order to evaluate the 
effect of increasing the horizontal resolution of the GLAS model 
and to assess the impact of sounding data v;ith this higher reso- 
lution model. Table 2 presents the results of changing from a 
4°x5° to a 2 l/2®x3® horizontal grid. It is clear from the table 
that the overall effect of increasing the horizontal resolution 
is to improve the quality of both searlevel and 500 mb prognoses. 

Table 3 presents the impact of satellite sounding data 
insertions on synoptic forecasts with the 2 1/2° x3° model. Com- 
parison of this table with Table 1 reveals that a more substantial 
beneficial impact from satellite data has resulted with the use of 
the 'higher resolution model. No major changes for the worse have 
resulted from the inclusion of satellite data while five sea-level 
and four 500 mb prognostic charts have been significantly improved. 
Once again the effect of satellite data is largest during the latter 
half of the forecast period. - . 

As an illustration of the impact of satellite-derived 
temperature soundings on prognostic flow patterns, the 72-hour 
forecast from February- 19, OOZ is included in this paper. This 
was a case in which large beneficial impact was also evident in 
the statistical scores. At the initial time of the forecast, 
the analysis showed a moderately intense low-pressure system 
located off the northwest coast of the United States. As this 
system moved inland, a new lov/ developed along an already existing 
stationary front and became the dominant feature. During the 
first 48 hours of the model forecast, only slight differences 
existed between the 2SAT and NOSAT sea-level prognoses as both 
correctly predicted the intensification and southeastward displace- 
ment of this low. After 48 hours, however, a radically different 
pattern emerged: the 2SAT forecast predicted the low to recurve 

and accelerate to the northeast, while the NOSAT forecast failed 
to move the low center. The NMC operational PE hemispheric fore- 
cast model also failed to move the low northeastward. This low 
was associated with both heavy snow and severe local storms during 
its displacement, and significantly different local weather fore- 
casts of those major events would have resulted from the different 
prognoses. 


Figures la and lb depict the 72-hour sea-level pressure 
prognoses from February 19 for the 2SAT and NOSAT cases, respectively, 
while the corresponding analysis is depicted in Figure Ic. Figure Id 
shows the NMC operational 72-hour forecast. A comparison of these 
charts reveals the very significant improvement in the displacement 
of the low-pressure system that resulted from the inclusion of 
satellite data. Substantial improvements in the displacement of 
the associated upper-level trough and vorticity maximum {charts 
not shown) were also observed. 
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Table 2. Results of a subjective assessment of NOSAT 

forecasts from the GLAS 2?2°x3'’ and 4'’x5® models. 


2a. Sea Level Pressure 


CATEGORY 

12h 

24h 

36h 

48h 

60h 

72h 

TOTAL* 

A 

0 

0 

' 1 

1 

0 

0 

2 

B 

2 

i 

6 

3 

2 

3 

4 

20 

C 

6 

3 

5 

6 

7 

3 

30 

D 

2 

1 

1 

1 

0 

1 

6 

E 

1 

0 

0 

0 

0 

0 

0 

0 


2b. 

500 mb Geopotential Height 



CATEGORY 

12h 

24h 

36h 

48h 

6 Oh 

72h 

TOTAL* 

A 

0 

0 

0 

1 

0 

1 

2 

B 

0 

2 

4 

4 

6 

1 

17 

C 

10 

7 

5 

5 

4 

5 

36 

D 

0 

1 ■ 

1 

0 

0 

1 

3 

E 

0 

0 

0 

0 

0 

0 

0 


* Total number of prognostic charts in each category. 


10 


























Table 3. Results of a subjective assessment of 2SAT and 
3ST0SAT forecasts using the GLAS 2%®x3° model. 


3a. Sea Level Pressure 


CATEGORY 

12h 

24h 

36h 

48h 

r - 

6 Oh 

72h** 

TOTAL* 

A 

0 

0 

0 

1 

2 

1 

5 

B 

0 

2 

2 

2 

1 

3 

10 

c 

11 

8 

7 

7 

5 

5 

• 43 

D 


' 1 

2 

1 

3 

1 

8 

E 

0 

>- 

0 

0 

0 

0 

0 

0 


3b. 500 mb Geopotential Height 


CATEGORY 


24h 

36h 

48h 

•60h 

72h** 

TOTAL* 

A 

0 

0 

0 

1 

1 

2 

4 

B 

0 

1 

3 

2 

3 

1 

.10 

C 

11 

9 

7 

7 

6 

5 

45 

D 

0 

1 

1 

1 

1 

0 

4 

E 

0 

0 

0 

0 

0 

0 

0 


* Total number of prognostic charts in each category. 

**3 forecasts missing. 
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a.) SAT Forecast 


b) NOSAT Forecast 


c) NMC Analysis for Verification 


d) IMMC Operational Forecast 


Sea level pressure plots for OOOOGMT 22 February 1976. Figure la shows 
the 72h forecast from SAT initial data at OOOOGMT 19 February 1976; 
Figure lb shows the corresponding NOSAT forecast; Figure Ic shows the 
verification analysis ■ (operational NMC analysis) ; Figure Id shows the 
NMC operational 72h forecast. Dots represent past positions of low 
at 12h intervals. 







Re ferences 


Atkins, M. J., and M. V, Jones, 1975: An experiment to determine 

the value of satellite infrared spectrometer (SIRS) data in 
numerical forecasting. Meteor, Mag ., 104 , pp. 125-142. 

Halem, M. , M. Ghil, R. Atlas, J. Susskind, and W. J. Quirk, 1978: 
The GISS Sounding Temperature Impact Test . NASA Tech. Memo. 
78063 , 4^ + xiii pp. 



QUANTITATIVE PRECIPITATION FORECASTS FROM THE GLAS MODEL 

(R. Atlas) 

The skill of the GLAS 4° latitude by 5° longitude 
general circulation model (GCM) model in predicting precipi- 
tation amount at U. S. gridpoints has been evaluated for 10 
of the forecasts that were generated during the February 
1976 DST period. Observed accumulations of precipitation 
at each of the 57 gridpoints located within the continental 
U. S. were obtained by first subjectively drawing isohyets 
on .facsimile charts of observed 24-hour accumulation at 
first-order U, S. stations, and then interpolating these 
observations to the model gridpoints. Since these observations 
cover the 24-hour period beginning at 1200 GMT while all of 
the forecasts in the GLAS sample were generated from 0000 GMT, 
only the 12-36 hour and 36-60 hour forecasts of precipitation 
amount were verified. 

Figures 1 through 4 present the model’s skill in 
forecasting measurable precipitation >,01 in., and precipitation 
amounts of >.l in. for each of the 24-hour periods. The 
three curves in the lower portion of each figure are: 

. PC - the number of gridpoints at which precipitation 
in a given category was forecast correctly, 

. PF - the total number of gridpoints at which precipitation 
in a given category was forecast, 

. PO - the total number of gridpoints at which preci-pi-tation 
in a given category was observed. 

The ratio PF/PO is the bias and represents the degree to which the 
model overforecasts or underforecasts the extent of precipitation. 
This is indicated by the distance between the PF and PO curves on 
the graph. FC/PO is called the prefigurence and measures the 
areal extent of precipitation that is correctly forecast in each 
category; FC/PF is the post agreement and measures the percentage 
of precipi’^ation forecasts which verify. These latter measures 
are represented by the distance between the FC and PO curves and 
FC and PF curves, respectively. 

The bias, prefigurence, and post agreement can be combined 
into a single ratio, known as the threat score (Th), and this 
latter quantity is plotted as the upper curve in Figures 1 through 4. 
The threat score is defined as follows: 

Th = FC/(PF + PO - FC) 

A threat score of zero represents no skill, while higher values 
denote increasing ‘skill. 
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FEBRUARY 1976 


Figure 1. Verification of 12-36 hour QPF for ^ .01". 
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NO. FORECAST 


NO. OBSERVED 


Figure 2. Verification of 36-60 hour QPF for ^ ,01". 


L6 





THREAT SCORE 

(AV. = .34) 

PSTAG = .44 

PRFIG = .59 
* 



Figure 3 . 


Verification of 12-36 hour QPF for > .1". 
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From Figures 1 and 2 , it can be seen that the GLAS model 
possesses a high degree of skill in forecasting the areas where 
measurable precipitation occurs (as measured by the- prefigurence) , 
but that a considerable bias toward overforecasting exists. The 
large degree of day-to-day variability of the threat score is due 
.almost entirely to changes in the amount of bias of the forecasts. 
The average threat score for the entire sample compares favorably 
with results previously reported by Druyan , et al . (1975). 

Figures 3 and 4 show that the GLAS model possesses 
considerably less skill in predicting precipitation amounts of 
>.l in. even though the bias of the forecasts has been substantially 
reduced. It is interesting to note that the loss of skill in going 
to higher precipitation amounts is greater than the loss of skill 
due to time, as can- be -seen from a^ comparison of the- threat scores 
in Figures 2 and 3. This loss of skill and decrease in overfore- 
casting continues as the precipitation amount category increases. 

The GLAS model failed to predict any of the gridpoints in this 
sample with precipitation amounts of >1 in. In this category PP, 

FC, and Th were equal to zero. 

All of the forecasts described above were generated 
from NOSAT initial conditions, in which no satellite-derived 
temperature soundings were included. When the forecasts were rerun 
from SAT initial conditions, with satellite soundings inserted, 
a small improvement in the average threat score of +3 percent ' 
occurred for the entire sample. A time stratification of the 
impact revealed a slight decrease in skill' of -2 percent in the 
12-36 hour predictions of >.01 in. and >.l in. combined, and a 
more significant increase in skill of +6 percent in the 36-60 hour 
predictions. This agrees with our. synoptic evaluations which 
indicated the maximum impact of sounding data over the U.S. to 
occur in the latter half of the forecast period. 

Reference 


Druyan, L. M., R. C. J. Somerville, and W, J. Quirk (1975): 

Extended range forecasts with the GISS model of the global 
atmosphere. Mon. Wea. Rev ., 103 , pp. 779-795. 
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DEVELOPMENT OP A COMPUTERIZED PROCEDURE FOR THE 
PREDICTION OF SEVERE LOCAL STORM ' POTENTIAL 

(R-. Atlas) 

The development of severe thunderstorms and tornadoes 
occurs on a much smaller scale than is currently resolved by the 
operational observational network or numerical weather prediction 
(NWP) models. However, these phenomena are usually imbedded within 
a larger region of instability and destabilization, and, therefore, 
the prediction of severe local storm potential depends to a large 
extent upon the ability to predict the synoptic scale setting 
favorable for intense convective development. 

If it is assumed that equivalent potential temperature 
(0e) is conservative following .the motion and let cr=90e/3z, 
evaluated between the 850 mb and 500 mb levels, represent con- 
vective instability (negative if unstable), then the following 
equation for local destabilization may be readily derived: 


(v-V0e) _ 3 (wg) 

3t 3z 

( 1 ) ( 2 ) 

This equation states that the local rate of change of convective 
instability depends upon the effects of differential equivalent 
potential temperature advection (term 1) and vertical stability 
advection (term 2). Accordingly, whenever an initially unstable 
environment is coupled with increasing upward vertical motion 
with elevation in the lower a’tmosphere, then -3L_(v-V0e) 

3z 

negative for convective instability to be increased or maintained. 
This occurs whenever there is a horizontal transport of warm, moist 
air at low levels of the atmosphere and/or cold, dry air at higher 
levels . 


Previous studies have shown the tendency for severe 
weather to develop in regions of strong differential advection and 
have suggested utilizing differential advection as a severe storms 
predictor. According to Whitney and Miller (1956), differential 
advection increases the potential for severe convective activity 
in three ways; { i) it decreases the amount of work required to lift 
low level air; -(ii) it decreases the height through which the air 
must be lifted before it will rise freely; and (iii) it increases 
the amount of energy that is available to a convective system. 

This process is often initiated by the movement of a cold trough 
toward the central United States (Newton, 1963) but may be present 
in a variety of synoptic situations. 
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The purpose of this investiqation is to determine the 
magnitudes of convective instability and differential advection 
associated with severe storm developments and to develop a 
computerized procedure for evaluating the potential for severe 
storms for a particular area. To this end, it is necessary to 
determine ( i) the percentage of severe local storm development 
in which differential advection plays a significant role, 

(ii) what parameters must be considered in conjunction with 
differential advection, (iii) with what resolution differential 
advection must be resolved, and (iv) how well current NWP models 
can predict differential advection. 

A preliminary version of a computerized severe storm 
model (CSSM) was developed at the outset of this investigation 
in order to provide a logical framework for evaluating the 
relative utility of various large-scale indicators of severe 
storm potential. The model uses a grid system of 4° lat. by 
5® long, to be compatible with the GLAS GCM, and covers the 
continental U.S. However, in this study only the 22 gridpoints 
east of lOS^W longitude were considered. This area was chosen 
for its high frequency of severe weather and to mitigate the 
effects of topography. A second version of the CSSM compatible 
v/ith the 2 1/2° latitude by 3° longitude version of the GLAS 
model is also being tested. In this case only the U. S. grid- 
points east of 102°W are considered. 

Input into the model consists of specific humidity, 
temperature, and wind velocity from either an objective analysis 
or numerical forecast tape. From these a storm potential index 
(SPI), differential advection (DA), convective instability (a) 
and twelve other measures of instability or destabilization are 
derived. The calculation of storm potential for this initial 
version of CSSM was based to a large extent upon a diagnostic 
study of 20 different severe storm outbreaks by Goldstein (1976). 

In this study it was shov/n-that in each of the outbreaks virtually 
all of the severe weather occurred within much larger regions 
where differential advection of -,5°C per hour between the 850 mb 
and 500 mb levels coincided with an initially unstable environment. 
Accordingly, the CSSM utilizes a and DA as its primary indicators 
of storm potential. However, it also attempts to delineate the 
severe weather areas more precisely by considering conventional 
indicators as well. It begins by testing the model atmosphere 
for convective instability and differential advection in combi- 
nation, and then tests for additional destabilizing effects due 
to vorticity advection and wind shear. On the basis of how 
these factors combine, a number of +1 to +12 for initially unstable 
situations or -1 to -12 for initially stable situations can be 
assigned to each grid box surrounding a model gridpoint, where 
±12 indicates the highest potential in each case. 

Two phases of experimentation have been partially 
completed. The first phase of testing of the CSSM has been to 
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conduct diagnostic studies to assess the potential usefulness 
of DA and a as specified by GLAS model, analyses. Since 
analyses c.an be cons-ideired as "perfect prognostic charts," this 
study can be used to determine the maximum predictive capability 
of the CSSM for medium-range forecasts (greater than 12 hours) 
or the actual predictive capability for short-range forecasts 
(0-12 hours) . 

Table 1 presents ( in abridged- form) some of the 
preliminary results of this study for the period February 2-29, 
1976. During this period 109 severe thunderstorms and tornadoes 
were reported within the CSSM’s region. From the table, it can 
be seen that negative differential advection occurred in a large ' 
percentage of the severe weather reports and that utilizing o 
and DA in combination is considerably more useful than convective 
instability alone in specifying areas of severe weather potential. 
Modifications to the original SPI calculation are currently being 
made to reflect these results. 


Table 1. Percentage of Storms Occurring 0-12 Hrs. 
After a and DA Analysis. 


DA 

- -10 

< -5 

1 A 
O 

< +5 

< +10 

TOTAL 

■1. 

21.1 

0 

11.0 

13.9 

9.1 

55.1 

-.5 

0 

4.6 

0 

0 

10.1 

14.7 

0 

0 

0 

1.8 

11.0 

1.8 

14.6 

+ .5 

0.9 

13.8’ 

0.9 

0 

0 

15.6 

TOTAL 

22.0 

18.4 

13.7 

24.9 

21.0 

100 


In the second phase of experimentation, the CSSM has 
been applied to actual forecasts from the GLAS GCM. As an 
example of the skill of this model in predicting a, DA, and 
absolute values of SPI (representing the combination DA with 
vorticity advection and wind shear) , the mean absolute errors 
from the 72-hour forecast from February 19, 1976, OOZ (using the 
2 1/2° X 3° version of the model) are presented in Figures 1 
through 3. This- period contained the largest of the severe 
storm outbreaks that occurred during February 1976 and also 
was a case in which a major improvement to the predicted develop- 
ment and displacement of synoptic scale features resulted from 
the insertion of satellite-derived temperature soundings in 
data sparse regions (Atlas, et al, 1978). 
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Figure 1. Mean absolute errors 
in the 72-hour prediction of 
convective instability (a) from 
February 19, 1976 OOZ. 
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Figure 2 . Mean absolute errors 
in the 72 -hour prediction of 
differential advection (DA) from 
February 19, 1976 OOZ. 


Figure 3. Mean absolute errors 
in the 72-hour prediction of the 
storm potential index (SPI) from 
February 19, 1976 OOZ. 









Curves N and S represent the forecast errors resulting 
from NOSAT and SAT initial conditions, respectively, while curve p 
represents a persistence forecast for the same time period. A 
comparison of -these curves in each of the figures reveals that 
both versions of the GLAS model possess considerable skill over 
persistence throughout the first 48 hours of the forecast, and 
that a significant reduction in the magnitude of the errors during 
the latter half of the forecast has resulted from the insertion 
of satellite-derived temperature soundings in data sparse regions. 
In this case the 48-60 hour predicted area of high severe storm 
potential generated from SAT initial conditions contained more 
than 70 percent of all reported severe thunderstorms and tornadoes 
while the corresponding NOSAT prediction accounted for only 
40 percent of the severe reports. 
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TROPICAL DISTURBANCES IN THE GLAS GCM SIMULATIONS 

(M. A. Es toque) 

An important aspect of general circulation modeling 
work that has not yet received sufficient attention is the 
verification of the model simulations of tropical weather 
disturbances. A 1-month study was conducted to examine the 
GLAS general circulation model (GCM) simulations of such 
disturbances. The area of study was limited to the tropical 
Pacific region, together with the so-called Indonesian Maritime 
Continent. This area extends from 90°E to 150°W longitude and 
from 25°S to 25®N latitude. The Indonesian region has been 
included because it is the strongest source of heat for the 
global atmosphere during the winter season. This heat source 
is, therefore, an important component for driving the general 
circulation of the atmosphere. 

Some of the interesting questions that could eventually 
be answered by a study such as that envisioned above arej 

What types of synoptic scale disturbances 
are generated by the GLAS GCM? What are their 
characteristic three-dimensional structures, wave- 
lengths, and speeds of propagation? What physical 
processes are responsible for their generation, 
maintenance, and dissipation? What large-scale 
flow parameters determine their wavelengths and 
phase speed? 

Partial answers to the above questions are available 
from previous observational as well as theoretical studies. 

These studies indicate that tropical weather disturbances may 
be classified primarily into three categories; wavelike 
disturbances, vortices, and linear disturbances. The waves are 
characterized, in general, by cold troughs in the lower half of 
the troposphere; they are most intense near the 700 mb level. 

The typical wavelength is about 3000 Km and the phase speed is 
about 7 m sec“^ westward. Cloudiness and rainfall occur most 
often to the east of the surface trough. Vortices may be of 
two kinds — low-level and high-level cyclones. Low-level cyclones 
are warm-cored in the entire troposphere and are most. intense 
in the middle and the upper troposphere. The best known example 
of the low-level cyclone is the hurricane. Upper level cyclones 
are cold-cored in the troposphere. They are most intense in the 
middle and the upper troposphere. In the last category, the 
linear disturbances are those whose cloudiness and rainfall 
patterns tend to be oriented along a line. A common example 
of a linear weather disturbance is the intertropical convergence 
zone (ITCZ). Another example is the shear line, which is simply 
an extension of a temperate latitude cold front into the tropics. 

The specific objective of the initial study was to 
make a general survey of the types of tropical disturbances that 
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occur in the GLAS GCM’ simulations. For this purpose, the 
results of a numerical integration for a winter season were used. 
Simulated winds at individual gridpoints were available at 12 
hourly intervals for various levels in the atmosphere. Other 
variables., such, as temperature., surface pressure, specific 
humidity, and vertical velocity, were also available. It was 
noticed that the gridpolnt distributions of these variables 
were characterized by large-amplitude , spurious waves with 
wavelengths equal to double the grid length. In order to 
eliminate these waves, Shuman's 9-point smoother was applied. 
Various types of analysis (e.g., streamline, time cross-sections) 
were done on the smoothed and unsmoothed array of gridpoint 
values. Some of the more interesting results of the analysis 
are presented below. 

The results indicated, in general, that the three 
categories of tropical disturbances described above occurred in 
the GLAS GCM simulations. The streamline fields at the 900 mb 
level often showed vortices and the ITCZ. An example of a 900 mb 
map is shown in Figure 1. One may see two vortices, a cyclonic 
vortex (longitude 135 °E, latitude 15 °S) and an anticyclonic 
vortex (longitude 150°E, latitude 4.*’S). It is interesting to 
note the unusual wind convergence towards the center of the 
anticyclonic vortex. In the same figure, one may see the ITCZ 
along an east-west direction just south of the equator over 
the western and eastern sections of the region. The orientation 
and the location of the ITCZ is rather realistic. Wave-type 
disturbances may be seen in Figure 2; this map represents 
conditions near the 700 mb level. The axis of the troughs are 
indicated by dashed lines. 

Although the results of the analysis presented above 
are very encouraging, there appear to be many unrealistic 
features of the simulations. One such feature is the abrupt 
changes in the wind with height. This is illustrated strikingly 
in Figure 3, which shows a vertical time cross-section over a 
gridpoint near the equator. In this diagram, the horizontal 
wind vector is drawn in the conventional meteorological sense, 
i.e., arrov7 pointing upward on the page means northward flow; 
downward -it-ans southward flow. Note especially the winds on 
the right half of the diagram. One can see changes from one 
direction to the opposite direction, going from Level 8 to 
Level 7. The same abrupt changes occur from Level 4 to Level 3. 
These anomalous changes in the wind with height may be due to 
numerical scheme for integrations. 

Another unrealistic feature of the simulations is con- 
cerned with the motion of disturbances. It has been noted that 
the simulated disturbances tend to be stationary or move slowly 
eastward; whereas, in the real atmosphere, the disturbances 
almost always move westward in the tropics. This unrealistic 
characteristic of the simulations is illustrated in Figures 4 
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Figure 3. Vertical-time cross-section of the wind. 










and 5. In these diagrams, the horizontal profiles {along 2®N 
latitude at the 300 mb level) of the northward wind component 
are plotted for successive days. Figure 4 represents unsnoothed 
profiles; Figure 5 shows the corresponding smoothed profiles. 

Note the overall lack of appreciable displacement of the wave 
forms along the east-west direction. There are other shortcomings 
of the simulations which cannot be described in detail in this 
brief report. One of the more glaring ones is the large 
variations (horizontal and time) in the surface pressure and 
the temperature. Another is the general lack of time and space 
continuity in the intensity and the notion of the simulated 
disturbances . 

In summary, the goal of the study was to make an 
initial determination of the realism of the GLAS GCM simulations 
of tropical disturbances. This was done by performing a synoptic 
analysis of the gridpoint values of various meteorological 
variables from a winter simulation. The analysis indicated 
that the GLAS GCM was able to generate various types of tropical 
disturbances — waves, vortices, and. linear disturbances — which 
exhibited many realistic features. On the other hand, the 
analysis also revealed various shortcomings in the simulations. 

It is desirable that this initial study be continued to make a 
thorough validation of the structure, phase speeds, and other 
characteristics of the simulated synoptic disturbances. The 
results of such a validation would be valuable in further 
improvement of the present GLAS general circulation model. 
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Figure 5. East-West profiles of the northward wind component (smoothed) 






TIME-CONTINUOUS ASSIMILATION OF REMOTE-SOUNDING DATA 
AND ITS EFFECT ON WEATHER FORECASTING 

(M. Ghil, M. Halem, and R. Atlas) 

Time-continuous methods for the four-dimensional assimilation 
of remote- sounding temperatures are presented in this paper to 
include time-continuous versions of direct insertion, successive 
correction, and statistical linear regression. They are applied 
to temperature sounding data obtained from radiance measurements 
performed by instruments aboard the polar-orbiting satellites 
NOAA-4 and Nimbus-6. The data were collected during the U. S. 

Data System Test in January-March 1976. 

A comprehensive series of experiments was performed to study 
the effects of using various amounts of satellite data and differ- 
ing methods of assimilation. The experiments included the assimila- 
tion of data from the NOAA-4 satellite only, from Nimbus-6 only, 
and combined data from both satellites. Other experiments involved 
variations in the application of the GLAS time-continuous statistical 
assimilation methods and of asynoptic successive correction 
methods. Intermittent assimilation was also experimented with, 
and its results compared with those of time-continuous assimilation. 

Atmospheric states determined in the assimilation experiments 
were then used as initial states for a number of 3-day forecasts 
corresponding to each experiment. The- effect of the satellite 
data was studied using the following criteria: (1) differences 
between the initial states produced with and without utilization 
of satellite data, (2) differences between numerical predictions 
made from these initial states, and (3) differences in local 
weather forecasts resulting from the large-scale numerical 
predictions. 

Initial state differences were evaluated in terms of magnitude 
and location of large-scale differences between meteorological 
fields. Numerical prediction differences were evaluated in terms 
of S-i skill scores and RMS errors, as well as by synoptic case 
studies. An automated forecasting model (AFM) based on quasi- 
geostrophic theory and on subjective forecasting principles was 
developed to facilitate the objective evaluation of differences 
produced in local weather forecasts, especially precipitation 
forecasts. 

Conclusions from these studies are that: (1) satellite- 

derived temperature data can -have a modest, but statistically 
significant positive impact on numerical weather prediction in 
the 2- to 3-day range; (2) this impact is highly sensitive to 
the quantity of data available and increases with data quantity; 
and (3) the assimilation method plays a major role in the magnitude 
of the impact obtained for the same data. 
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RECURSIVE ESTIMATION FOR NUMERICAL WEATHER PREDICTION 


(M. Ghil and J. Tavantzi-s-) 

The initialization problem of numerical weather pre- 
diction (NWP) consists in providing the initial data necessary 
for solving the -prognostic systems used in NWP, These data are 
not provided with the required completeness or with sufficient 
accuracy by existing observing systems. Hence, the need to 
compute the data at initial time, t = 0, from incomplete and 
error-contaminated observations at preceding times -T < t < 0. 
Such observations are provided by a variety of unconventional, 
asynoptic observing systems, in particular by polar-orbiting 
and by geostationary satellites; these observations add informa- 
tion to that obtained by the conventional, ground-based synoptic 
meteorological network. Given both asynoptic' and synoptic data, 
GLAS proposes to obtain optimal estimates of initial states for 
numerical forecasts using a systematic approach to four-dimen- 
sional (4-D) data assimilation. The approach is based on the 
control-theoretical methods of filtering and of optimal estima-. 
tion. These methods, commonly associated with the work of 
Kalman, Bucy and many others, have proven remarkably successful 
in a number of engineering problems. Attempts have also been 
made to apply, them to problems in dynamic meteorology by Miyakoda 
and Talagrand (1971), Petersen (1968, 1970), Phillips (1971), 
and Talagrand and Miyakoda (1971). 

The application of recursive estimation methods to a 
simple model exhibiting important features of large-scale 
atmospheric flow will be studied-, namely to a shallow-fluid 
model. In a rotating Cartesian x ,y-coordinate system, the 
equations governing the model are 
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4)y + fU = 0, 
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u4>x 
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+ 

d)(Ujj+Vy)= 0; 


( 1 ) 


here u and v are the velocities in the x and y directions, 
respectively, and <|) is the geopotential (Ghil et al . , 1977). 

The domain of definition is a rectangle in the x,y-plane 
with u, V, cj> periodic in x , so that the domain represents 
in fact a zonal channel. 

t 

For system (1) the estimation problem can be formulated 
as follows; the state of the system' w = ((j),u,v) is approximated 
at t = -T and the statistical structure of the error in this ■ 
approximation is assumed to be known. Then some of the variables 
are observed continuously or discretely with a certain error 
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in the interval -T < t < 0 . The best estimate w = ($,u,v) 
of w at t = 0 is sought. This estimate will then be used as 
the initial state for a numerical forecast starting at t = 0. 

The effect of different observational configurations 
and of various error structures on the results of the optimal 
estimation will be studied. . The approach will then be compared 
with previous methods for 4-D data assimilation: updating 

(Charney et al . , 1969) , ^optimal interpolation (Rutherford, 1972; 
Ghil et al . , 1979a) and others. 

The first model considered is a linearized g-plane 
approximation of ( 1.) ; 


u^ + UUjj + - fv = 0 

v^ + UVjj + fu = 0 (2) 

<1)^ .+ _fUv_ = 0_ 


Assuming u, v and are independent of y, the domain is the line 
interval which represents the 45 °N latitude circle, and requires 
periodicity at- the end points.' Therefore w=(u,v,cj)). For 
geostrophic initial conditions, i.e.. 


w(x,-T) = (0,d)jj/f,<^) , (3) 

the solution w{x,t) to (2,3) contains only slow waves. The 
observed initial condition at t = -T is (3), on which normally 
distributed noise is superimposed, is assumed. The solution 
w(x,t) is further observed during the assimilation period 
-T < t < 0 at discrete values of x and t. These observations 
are also made in the presence of noise. GLAS is interested in 
estimating w(x,t) at t = 0 , given that new information is made 
available; in particular, GLAS is interested in estimating the 
slow waves. 


A Richtmeyer two-step scheme is used to solve (2) 
numerically and express the discretized equations in the form 


^n+1 = 'i'n ^n ' • <4) 

where Wj^ is the vector obtained by evaluating w(x,t) at the 
grid points iAx at time -T+nAt, and is the numerical 
amplification matrix. 
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Let be the vector of observations at 

= -T+(n+l)At7 It is expressed as 


' Sn+l ^n+1 ~n+l ^n+1 ' 

where v bhe observational noise. The observation matrix 

reflects the fact that in general it is not the state varia- 
bles Wj^ themselves which are observed, and certainly not all 
of them; thus the^dimension of the vector is smaller than 
that of w . Let be the estimate at the (n+1)®^ step, and 

^n+1 covariance error matrix, 

Pn+1 = E(«"n+1 “ ?n+l>^^n+l " 5n+l>^) ' ' (6) 

with E the expectation operator and ( the matrix transpose 

operator. 

The problem is now to estimate Wj^+2 minimize 

the trace of P^+l^ i.e., to minimize the expected least-square 
error in the estimate. The solution to this problem is given 
by the Kalman filter. 

^ After an observation at time nAt, the updated estimate 

Wn(+) of the solution is given by a certain linear combination 
or the preliminary estimate and of the observation z^. 

This combination is chosen so that the estimation is (a) unbi'ased 
and (b) minimizes the trace of the covariance error matrix (6) 
(Kalman and Bucy, 1961) . These two conditions determine 
the Kalman Gain Matrix at the n“^“ step (Gelb, 1974)' 


.here is the error covariance matrix, the observation 
matrix (Eq. (5)) and the covariance noise matrix, Rj^=E(Vj^ 

The error covariance matrix update is given by 


?„( + ) = (I-Kj^Hn)Pn{“) 


(8a) 


The state estimate update now becomes 



The error covariance matrix Pj^ evolves in time according to 


Pn+1 (-) = 'i'n Pn ( + ) '^'n 


T 


( 10 ) 


In the absence of observations at time tj^^ = -T+mAt, 
the updates of the state and of are given by 


Pin ( + ) = Pm 


(3b) 


and 

Wj^ ( + ) = Wj^ (-) , ' (915) 

with 

(-) = Vl "m-1 ( + ) • (4") 

The Kalman filter technique was applied to the advection 

equation 

u^ + Ujj = 0 , 0<x<10 , 0<t (]^ 1 ) 

« 

with a triangular wave as "correct" initial data 

u(x,0) f(x) = X . , 0 < X < 5 , (12a) 

=10-x , 5<x<10 , (12b) 


and periodic boundary conditions 


u(0,t) = u(10,t), 0 < t 


(13) 


The initial data were given with a uniformly distributed 
random error of amplitude 0.1 Uj_^ , Later observations 
were made with a random noise or 0.01 u^„„.. 

One numerical computation in which the filter was 
applied to this situation is illustrated in Figures 1 and 2. 

The interval 0 < x < 10 has been partitioned into 10 subintervals 
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X£ = iAx , Ax = 1; the time step is At-= 1. The numerical 
scheme was upwind differencing (first order). The three arrows 
in Figure 1 indicate the location of observation points, 
x° = (2,5,8-).; i-n th-i-s example, observations were made at every 
time step. The situation thus corresponds to synoptic observa- 
tions being provided to a quasigeostrophic model with large 
time steps by a sparse conventional network x° . 

The. initial error is given by curve O* The abscissa in 
Figure 1 is x and the ordinate represents the error u - u between 
the true value u of the solution to (11-13) and its estimated 
value u. One can easily see that, after three time steps, the 
error becomes negligible. It is clear that the error u - G is 
reduced as the triangular wave f(x-t) travels through the 
"observational network" x®. 

Figure 2 depicts the behavior of the relative L 2 -error 
(RMS), I I u-u I 1 2/1 I u I I 2 f and of the weighted L 2 ~error given by the 
trace of the error covariance matrix, tr P. Both error measures 
decrease rapidly and reach a quasi-steady value after three time 
steps. The technique worked well in this simple case, yielding 
estimates of the solution u with errors not exceeding those in 
the. observations. 

The Kalman filter has been programmed for the system 
(2,3). The program is structured as a package to allow the 
investigation of Kalman filtering for different observation 
patterns. One test has been carried out in which the computa- 
tional network included 16 equidistant points. The .observational 
network was x°= (x^ ,Xo ,x^ ,x^ ,xg ,x^q ,x^ 2 ) w=(<|),u,v) was 

observed at each x^ belonging to x° every third time .step t^°=3nAt. 
This situation is similar to the one reported on for the advection 
equation (11-13). In this case rapid error decrease occurred 
also, so that after two observations, or 6nAt, the error fell 
from its initial high value of about 10 percent to the low value 
of observational error of about 1 percent. 

Preliminary results indicate that recursive estimation 
using the Kalman filter is very effective in 4-D data assimila- 
tion for the linearized shallow-fluid equations. GLAS intends 
to study observation patterns that simulate closely both the 
synoptic and the asynoptic observational network in the mid- 
latitudes of the Northern Hemisphere (Phillips, 1976). The 
present technique will then be compared to a number of pre- 
viously implemented assimilation methods. 

It will be especially interesting to compare recursive 
estimation to the linear regression methods which, under the 
name of "optimal interpolation," have recently attracted con- 
siderable interest in the meteorological community (Bengtsson, 

1975; Rutherford, 1972). The latter methods are purely 
statistical in nature, ignoring the dynamics of the atmosphere, 
as expressed in the equations of motion used in NWP. A step 
in the direction of combining dynamics with statistics appears 
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Figure 1. Evolution of pointwise estimation error u-u. 

Arrows^indicate observation points, x=2,5,8. The 
curve(^ indicates u-u after the n“?^h observation, 
u(iAx/nAt) - u(iAx, nAt) . 
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Figure 2. Evolution of two measures of estimation error; the 

relative L 2 “error | | u-u| | 9 / ||ull 2 and the weighted 
L 2 “error, trP = (Uj^-Uj^) 2 . 


39 


in the time-continuous statistical assimilation method (SAM) 
of Ghil etal. (1979a, b). 


GLAS shall attempt to -devise improvements in practical 
assimilation methods by an increased understanding of recursive 
estimation for simple prognostic systems. As a first step 
towards that goal, GLAS plans to extend the results to the non- 
linear system (1) by applying a generalized Kalman filter 
(Seinfeld et al . , 1970)'. 

Control theory should provide also a unified framework 
for both data assimilation and initialization methods. Filtering 
of high-frequency inertia-gravity noise can be performed within 
the same procedure, which assimilates asynoptic and synoptic 
observations (Ghil, 1979). This could eliminate the current 
practice -of following an assimilation cycle by a computationally 
costly initialization- calculation. 
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PARAMETERIZATION OF THE FRICTIONAL TURNING ANGLE 
FOR THE HOFFERT-SUD PLANETARY BOUNDAR-Y LAYER' MODEL 

(M. Helfand) 

The Hoffert-Sud Planetary Boundary Layer Model ( J. Atmos. 
Sci., 1976) is a similarity model in which the vertical coordinate 
is normalized by the boundary layer depth, h. The structures 
predicted by the model can be parameterized in terms of the non- 
dimensional variables, Rij^ = h/L, h/z^, and h^^/h, where L is the 
Mon in- Obukhov length, Zq the surface roughness length, and h =.35 u/f 
is the Ekman depth. The model predicts a drag coefficient, Cj^, a 
heat transfer coefficient, Cjj, and a frictional turning angle, 'Fp, 
as a function of these three similarity parameters. In the past, 
only Cjj and have been used in the parameterization implemented 
in the GLAS GCM. However, it is necessary to compute the turning 
angle in order for the GLAS GCM to assimilate sea surface wind 
observations, such as those soon to be available from the Seasat-A 
scatterometer . 

In order to influence the GCM forecast, the level 9 
wind field must be altered, to reflect the surface wind information. 
This requires that the turning of the wind between the surface and 
950 mb be accounted for. 

Figure 1 shows i'p as a function of Rij^ for four values 
of Z^/h and two values of n^/h as determined by numerous inte- 
grations of the Hoffert-Sud model. An analytic approximation to 
this function has been constructed: 'Fp is taken as the log linear 

function of Rij^ for R-ij^ <-3 and as an exponential function of Rij^ 
for -3<Rij^<0. For the stable case (Rij^ >0) the Richardson number 
is set equal to zero. Unfortunately, integrations of the Hoffert-Sud 
model imply that Yp increases for decreasing h /h for the neutral 
case, behavior that is not physically reasonable. The reason for 
this anomalous behavior at neutral stability is currently under 
investigation. 
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TURNING ANGLE ijf 



Figure 1. The frictional turning angle Xf a function of bulk Richardson 
nuinber for several values of h^/h and z^/h. The solid lines 

represent curves for hc/h = 1/ while the dashed lines represent 
curves for h^,/h = 2. The value of z^/h is presented next to 


FURTHER DEVELOPMENTS OF A PLANETARY BOUNDARY LAYER PARAMETERIZATION 


(D. Randall) 

A reailsfic' parameterization of the planetary boundary 
layer (PBL) must have the following features: 

A. An accurate method for determining the time 
evolution of the PBL depth during periods of 
both deepening and shallowing; 

B. A method of calculating the turbulent fluxes 
throughout the depth of the PBL/ under all 
conditions of stability, cloudiness, etc,; 

C. A comprehensive parameterization of the 
interaction of the PBL with clouds; 

D. A simple and complete method of incorporating 
the parameterized PBL into the framework of a 
numerical prediction model. 

Each of these major points was addressed by Randall (1976) in a 
parameterization for the six- and twelve-level UCLA general 
circulation models. Development of the parameterization is 
continuing, however, and all of these features of the para- 
meterization have been or are being significantly revised. The 
purpose of this report is to summarize the current status of the 
parameterization. 

The Entrainment Parameterization 


Because the PBL is defined as the layer within which 
the turbulence is a dominant process, it is natural to discuss 
the deepening and shallowing of the PBL in terms of the conservation 
of turbulence kinetic energy (TKE) . In vertically- integrated form, 
this conservation law is approximately given by 


^^M B + S Bj^5Pj^/’g, ( 1 ) 

where E is the entrainment mass flux; ej^ is the vertically-averaged 
TKE density; B and S are the rates of buoyant and mechanical pro- 
duction of TKE, respectively; is the vertically-averaged 
dissipation rate; <Spj^ is the PBL depth in terms of pressure; and 
g is gravity. The gravity wave radiation and cumulus production 
of TKE are neglected; and it is assumed that the time rate of change 
.of the TKE content, ej^6pj^, is dominated by the time rate of change 
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of 5pJ^^ . Por convenience, a turbulence velocity scale a is 
defined by 

where pj^ is the vertically-averaged density. 

It is assumed that 

= aic^, (3) 

and that 

pMcT^ "" ^2^' ■ (4) 

where P is the gross rate of TKE production, i.e., the integral 
of the production rate over those levels where it is positive. 

The constants a, and a 2 are the only empirical parameters of the 
theory. The values aj = 0.75 and 32 = 0.84 are tentatively assigned. 

The simple theory sketched above takes into account the 
production of TKE at all levels within the PBL, including possible 
cloud layers. It also allows shallowing of the PBL, by permitting 
negative entrainment whenever the dissipation rate exceeds the 
net production rate. 

The en’trainment theory has been tested using a simple 
mixed layer model. The results (Figure 1) show a realistic 
diurnal variation of the PBL depth over land. 

Revision of the Turbulent Flux Parameterization 

Randall (1976) did not propose a new parameterization 
of the turbulent fluxes, but simply adapted the method of Deardorff 
(1972). Yamada (1976), Arya (1977), and others have proposed some- 
what more sophisticated theories, which are nevertheless generally 
similar to Deardorff s. But none of these theories takes into 
account the effects of stratocumulus or cumulus clouds, which may 
drastically enhance the turbulent fluxes by dominating the pro- 
duction of TKE, This is a particularly serious problem for strato- 
cumulus convection, because it often occurs over cool oceans where 
conventional theories incorrectly predict very weak turbulent 
transfer. The existing theories also behave badly in certain 
special but important limiting cases, such as free convection or 
equatorial flow. Also, although these theories predict the 
existence of characteristic interior structures for the PBL in 
various parameter ranges, and- so imply the existence of correspond- 
ing characteristic interior flux profiles, the nature of these pro- 
files has not yet been adequately, documented . Finally, all of these 
theories are for the equilibrium structure of the PBL, and the 
important problem of how near-equilibrium is maintained against • 
the inevitable disturbing influences has not been solved. 
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A new theory which overcomes these difficulties is 
currently being tested. For the purpose of determining the 
similarity functions whose existence is predicted by the theory, 
a "higher-order closure" model has been developed and thoroughly 
tested by Joseph Scire, a graduate student at MIT. After the 
similarity functions have been determined, the next step is to 
test the theory off-line. Only then will it be incorporated into 
a GCM. 

A Parameterization of the Layer Cloud Instability 

As discussed by Randall (1976), under some circumstances 
a stratocumulus cloud deck is susceptible to penetration by con- 
vective downdrafts whose nega-tive buoyancy is maintained by the 
evaporation of cloud' liquid water. Recently, it was discovered 
that the effects of this instability on the cloud layer can be 
parameterized on the basis of the assumption that the downdrafts 
reduce the fractional cloudiness in such a way that the remaining 
partly cloudy layer is neutrally stable with respect to the 
downdrafts. The theory not only yields the turbulent flux profiles 
in the unstable cloud layer, but also, as a bonus, provides a 
measure of the fractional cloudiness. The latter is predicted to 
decrease in a natural way as the degree of instability increases. 
This simple parameterization has already been tested off-line and 
programmed into the UCLA GCM. 

Testing of a Potential Enstrophy Conserving Scheme for Use With the 
the Parameterization 


In order to incorporate the PEL parameterization into a 
GCM, a generalized a-coordinate. system (Figure 2) was adopted in 
which both the earth's surface and the PEL top are coordinate 
surfaces. The PEL then consists of one or more model layers. 
Although the height and pressure of the PEL top vary in space and 
time according to a prognostic equation, the value of a at the 
PEL top is always the same. 

The quasigeostrophic dynamics of large-scale atmospheric 
flows are governed by the conservation of potential vorticity. For 
the PEL, the potential vorticity is inversely proportional to the 
PEL pressure depth, which can change very rapidly from one gridpoint 
to the next. Therefore, special care is needed to avoid large 
truncation errors in the advection of PEL potential vorticity 
if the quasigeostrophic dynamics are to be correctly simulated 
for the lowest model layers. 

Fortunately, a suitable scheme has recently been 
developed by Arakawa and Lamb. Tests of the scheme, together with 
the PEL parameterization, are currently in progress at UCLA. Eut 
it is desirable to have a "clean" test of the scheme as well, i.e., 
a test in which the new scheme is the only new feature in the 
model. With the assistance of Dr. Kingtse Mo, the GLAS GCM has 
been modified to provide such a clean test, and debugging is now 
underway. First results should be available in the near future. 
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Figure 2. The generalized u-coordinate . 
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The GLAS goal is to incorporate the. features outlined 
above into a comprehensive PEL parameterization for GCMs. 

GLAS is working toward this goal by systematically testing 
each feature separately. The individual features will be 
incorporated into the GLAS GCM as each becomes ready. 
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A FOURTH-ORDER FORECASTING MODEL 


(E. Kalnay-Rivas, D. Hoitsma, and P, Anol-ick) 

The GISS fourth-order model {Kalnay-Rivas, et al . , 

1977), which is a fourth-order, energy-conserving GCM on an 
unstaggered grid, had shown promising capabilities. It produced 
forecasts that showed an improvement over the second-order GISS 
forecasts with the same fine grid (4° x 5°) resolution, but that 
were somewhat inferior to the "ultrafine" forecasts. However, the 
first version of the model required excessive amounts of computer 
memory and time for execution. 

The model has been reprogrammed into the "fourth-order 
band model." The new program solves the primitive equations one 
latitude band at a time. The arrays are stored in an interlaced 
way, with all arrays being updated at the same latitude stored 
contiguously, and similarly for all arrays used in the computation 
of the time derivatives. This design of the program makes effective 
use of the virtual memory capability of the Modeling and Simulation 
Facility's IBM 370/165 or Amdahl computers. The virtual memory 
facility permits the execution of programs whose core size is larger 
than the one available, by placing the excess on disk and reading in 
those pages of information needed in the current calculations. The 
band fourth-order data structure and computations were constructed 
to optimize this virtual I/O process; a possible improvement may be 
to interlace also the arrays being updated with those used to compute 
the time derivatives. The use of the virtual memory facility avoids 
the explicit I/O used in the current Kern model, and yields a simpler 
program. (The band structure was suggested by G. Russell#) 

The band fourth-order model computations' have been 
optimized so that each time step is computed in the Amdahl in half 
the time required by the old model. The array structure has also 
been designed to reduce the amount of overall storage and high- 
speed memory by a factor of two (see Table 1). 


Table 1, Comparison of Fourth-Order Model Computing Requirements, 



Original 

Band 


4th-Order Model 

4th-Order Model 

Core (bytes) 

3500K 

1500K 

CPU time per step 



.( COMF1, COMP 2) 

34 sec 

17 sec 
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The model has been programmed so that it can use both 
the Matsuno and the leapfrog time schemes, the latter with the 
Robert time smoother. If and D(O^) denote the fields at time 
nAt and the corresponding time derivative computed from the space 
differences, the smoothed leapfrog scheme is 

Qn+1 ^ qH-I + 2 At D(Q^) 

+ ,5v(q”"1'+ - 2Q^) ■ 

with Q° = O®. The use of a smoothing coefficient requires a 
slightly smaller time step. For example, with v = .1, the 
model is marginally stable with At =- 288 sec.^ -compared to • 

At = 300 sec. for the Matsuno and leapfrog schemes. 

The smoothed leapfrog scheme is further modified to 
include the subgrid "physics" terms, and scaling and spatial 
smoothing procedures. Since the "physics" is called every 
few time steps, unless the leapfrog scheme is restarted after 
every call to the "physics," only one of the two consecutive 
fields will be affected. The restarting procedure is time- 
consuming, so it has been replaced by the following algorithm 


1 — 1 

1 

c 

= 

- 



-- +' 

S 

n 

n-1 

= qH-I 

+ 



where corresponds to the "physics" terms. This procedure, 
which ensures that the physics is applied to two consecutive 
time steps, .has been tested with good results. 

Preliminary Results . The first numerical integrations performe'd 
with the new band fourth-order model show dramatic improvements 
over forecasts made with the GEAS model with the same resolution 
("fine grid" or 4® latitude by 5® longitude). The quality of the 
forecasts is now comparable with those produced with the "ultrafine" 
(2.5® latitude by 3® longitude) version of the model. Results 
from a 3-day numerical integration are presented in Figure 1. 

During an extended 8-day integration of the new fourth-order 
model, the atmospheric systems remained remarkably smooth, 
exhibiting a realistic behavior both with respect to position 
and intensities. 
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Figure la. Observed sea level pressure (February 22 , 1976 - OZ) . 
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Figure lb. 72-hour forecast - second-order GLAS model. 






Figure Ic. 72-hour forecast - fourth-order GLAS model. 




THE EFFECT OF TIME SCHEMES ON MODEL SIMULATIONS 


(E. Kalnay-Rivas) 

Numerical experiments with simplified models .indicate 
that the choice of a time scheme can have an important effect 
on the model limit of predictability. For example, 20-day 
integrations with a quasigeostrophic model with parameterized 
radiation, friction and orography were carried out using the 
Matsuno and Lorenz four-cycle schemes, both with time steps 
of 1 hour and 1.5 hours. Numerical integrations starting from 
the same initial conditions and using different time schemes 
showed significant differences after 2 to 4 d'ays . After 20 days 
the only similarity between the two forecasts was in the long 
stationary wave forced by orography. There were significant 
differences even in the "climate" obtained through a 20-day 
average. It was found that th’e differences between forecasts 
run v;ith the same scheme and different time steps v;ere much less 
significant than those between different time schemes. 

In order to study the properties of efficient time 
schemes a series of experiments was performed with a simple 
shallow water equations (SWE) model. The Matsuno (M) scheme, 
which is currently used in the GLAS model, was modified by 
using updated pressures in the pressure gradient term of the 
momentum equations. This was done either in the corrector 
step (modified Matsuno I or MI), in the predictor step (Mil) 
and in both steps (Mill). The properties of modified leap 
frog schemes were studied, in which the pressure term 
computed at time nAt, is replaced by fi) (•(j)’^'*'^ ■ + /2 

(modified leapfrog I or LPI); (ii) (LFII ) ; and 

(iii) (<|)^“^ + 2(j)*^ + (LFIII). This last version 

was proposed by Shuman. Linear analysis indicates that if the 
time step is chosen for maximum stability, as is the case in 
practice, then HI and Mil can use a time step /r- larger than 
the regular Matsuno scheme, and they have better stability 
properties at high frequencies. Mill has worse stability 
properties than MI or Mil., Linear analysis also indicates that 
the TASU scheme (time alternating, space uncentered), formerly 
used in combination with the Matsuno scheme in the GLAS model 
(Arakawa, 1972) is not spacially symmetric and introduces 
numerical horizontal diffusion with a coefficient c^At , 
half as much as the Lax-Wendroff scheme. 2 

Some of the results of the SWE model are the following 
( i) MI and Mil can be used with a time step about 50 percent 
larger than the regular Matsuno scheme; (ii) of the modified 
leapfrog schemes, the only stable one is LFIII (proposed by 
Shuman) , which ran with At about 50 percent larger than the 
regular leapfrog. It can be combined with time smoothing 
(Robert filter) only for very small smoothing coefficients. 

The combination of LFIII with a small amount of time smoothing 
and the periodic application of a Shapiro filter produced 
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excellent results at a rate three times faster than the 
regular Matsuno scheme; (iii) Lorenz's N-cycle scheme, v/hich 
has very attractive properties, and is most stable with IS! = 4, 
is not much more efficient than the Matsuno scheme because it 
requires smaller time steps for stability. In agreement with 
linear analysis, the size of the maximum time step is very 
sensitive to diffusion terms. 

The modified Matsuno I (MI) scheme has been introduced 
into the GLAS model, where it has allowed an increase in the 
time step from 10 min. to 14 min. A 3-day forecast with the 
new scheme and At = 14 min. was virtually identical to the 
one obtained with the regular Matsuno scheme' and' At = 10 min. 

Semi-implicit schemes for the shallow water equations 
can be easily designed and solved (Kwizak and Robert, 1971) . 

For the baroclinic primitive equations, since there are as 
many inertia-gravity wave modes as vertical levels, semi- implicit 
schemes become rather complex and require considerable extra 
storage, which may offset the advantages of a larger time step. 

A semi- implicit scheme was designed in which only the external 
mode is written implicitly to take advantage of the fact that 
the phase speed of the Lamb wave , or external inertia-gravity 
wave mode, is c=±/YRT --300 m sec“^, and about 10 0 m sec~^ 
or less for internal modes. In a-coordinates this becomes a 
simple scheme for the following reason: the effect of the 

external mode enters the momentum equations mainly through the 
■a^ViT pressure term, where ir is the surface minus the top pressure. 
Linear analysis indicates that this term alone accounts for 
c^= ±/RT, or 85 percent of the total speed of the external 
mode. Therefore, a simplified semi- implicit scheme can be defined 
in which only the Vir term in the momentum equations, and the 
/q^V (m^)da term in the continuity equation are written 
implicitly. This scheme is essentially as simple as the 
corresponding one for the shallow water equations: it requires 

only the solution of a two-dimensional Helmholtz equation, and 
is uncoupled from the temperature equation . Simplified semi- 
implicit schemes were designed for Arakawa's schemes B and C 
equations and for the fourth-order GISS model on an unstaggered 
grid {Kalnay-Rivas , et al . , 1977). 

In preparation for the use of the simplified semi- 
implicit schemes inGCMs, several experiments were performed 
with the shallow water equations model. The semi- implicit scheme 
was combined with the splitting method used by Sundstrora (1974), 
which allows the reduction of the Helmholtz equation to a set of 
one-dimensional equations while still retaining formally second- 
order time truncation errors. 

The experiments indicate that whereas the regular 
leapfrog scheme was marginally stable with a time step of 15 m-i-n-.--. 
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and the LFIII was stable with a 20-minute time step, the semi- 
implicit schemes could use much larger steps, up to 50 minutes. 
However, the results using the split semi- implicit scheme showed 
serious distortions after a few days, even with time steps as 
short as 20 minutes. When successive overrelaxation was used, 
the results showed no distortions until after about 8 days, 
even with large time steps. An iterative scheme developed 
by P. Frederickson, in which the norm of the error is reduced 
after each iteration by a factor that is independent of the 
number of gridpoints or their closeness, is being tested for 
the solution of the Helmholtz equation. It should be noted 
that no checkerboarding pattern was observed, even though 
four .sets of Helmholtz equations were solved each time step, 
probably because the momentum equations in flux form were used, 
thus providing a weak coupling between adjacent points. 
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MODEL-SIMULATED TRANS IE^7T EDDIES IN THE NORTHERN HEMISPHERE 

(S. Srivatsangan) 

An -extensive study of non-stationary eddies in the 
GLaS nodel output is being conducted by the GLAS Modeling and 
Simulation Facility. This study is parallel to Blackmon (1976) 
and Blackmon, et al . , (1977), which shall hereafter be cited 
as Bl and B2, respectively. B1 and B2 are based on a 10-winter 
set of meteorological observations. 

In Bl and B2 the time series for such variables as sea 
level pressure and 500 mb geopotential height are first converted 
to time series of transient eddies of these variables by removing 
the time averages. Then the seasonal cycle is removed from the . 
transient-eddy data by subtracting either the first four harmonics 
of the annual cycle or the parabolic fit of each winter's data. 
Finally, the eddies are separated into three types, which are: 

( i) low-pass filtered or slow-transient eddies of periods from 
10 to 90 days; (ii) band-pass filtered or fast transient eddies of 
periods from 2.5 to 6 days; and (iii) high-pass filtered fluctuations 
with periods of less than 2 days. 

In this study the filter coefficients reported in Bl were 
used. The first four harmonics of the annual cycle were considered 
to be the seasonal cycle. The reader may note that the data are 
for single seasons only,' while the analyses in Rl and B2 are for 
10 seasons. 

The model integrations studied are called D091 and DlOO, 
both being for winter. D091 is characterized by a jump in sea 
surface temperature from month to month, whereas DlOO has a 
linearly interpolated sea surface temperature from day to day. 

DlOO integration was also performed using a filter for the two-grid 
wavelength (or 2Ax) eddies in the zonal direction, effective 
from approximately 60 °S to 60 °N. In other respects, dynamically, 

DO 91 and DlOO are similar to earlier GLAS model integrations. 

The results of model integrations are available at the 
o levels of _..e model. These are converted to pressure-level 
values and a time series data set for further computations. 

The results are in two parts. The first part contains 
the statistics for the 500 mb geopotential height, which can be 
compared with Bl , and the second part is for the winds at 
500 mb, which are compared with B2 . 

Figure 1 contains the temporal root mean square (RMS) 
values of the unfiltered 500 mb geopotential height eddies, both 
for model simulation D091 (January and February) and for atmospheric 
observations (Bl). The features in mid-latitudes are well defined 
in the model, and the maxima over the oceans are comparable in 
magnitude to atmospheric observations. Only the maximum over the 
U.S.S.R. is of considerably smaller value than in the observations. 
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(la) 


(lb) 


Figure 1. Temporal RMS values of 500 mb geopotential height 
eddies for GLAS model integration D019 (la) >and 
from Blackmon (lb) (1976). Units: geopotential 
meters . 



(2a) 


(2b) 


Figure 2. Temporal RMS values of 500 mb geopotential height 
eddies — low-pass filtered component — for GLAS 
model integration D091 (2a) and from Blackmon (2b) 
(1976). Units: geopotential meters. 
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Figure 2 has. the RMS values of the low-pass filtered 
500 mb geopotential height eddies. Here the maxima in the model 
are of smaller value in relation to the observations. The maximum 
over the U .-S-.-S-. R . / as sh'oWn in the observations,, appears only as 
an increase in value in the model. 

Unlike the RMS low-pass filtered values, the band-pass 
RMS values, given in Figure 3, are over the Pacific and are larger 
in the model than in the observations. The maximum over the 
Atlantic is comparable to the observed. The orientation of con- 
tours is also in consonance with atmospheric observation. 

^~rhe rest of the study's results deal with the statistics 
for the wind at 500 mb, as obtained in the integration DlOO for 
January and February. 

The seasonal mean zonal wind is represented in Figure 4. 
The jet stream maxima' straddling the east coasts of North America 
and Asia and the maximum over the Sahara and Arabia are almost 
identical to observations. In the polar latitudes a wavenumber 3 
pattern prevails, as in B2. 

One of the striking results of B2 is the location of the 
band-pass RMS maxima downstream and poleward of the Asian and North 
American jet streams. Figure 5 shows the band-pass RMS values 
of the meridional wind at 500 mb. The maxima in Figures 4 and 5 
show clearly that this is also observed in the GLAS model. 
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Figure 5. Temporal RMS values of band-pass filtered 
500 mb meridional wind for GLAS model 
integration DlOO (5a) and from Blackmon 
et al . (1977) (5b)* Units: ms“^. 



ATMOSPHERIC ABSORPTION SPECTRA NEAR 2200 CM ^ and 2400 CM ^ 


(T. Ho and J. Susskind) 

The spectral regions near 2200 cm“^ and 2400 cm“^ , 
on either side of the 4.3ym CO 2 band, are well suited for 
temperature sounding of the lower troposphere. Quantitative 
interpretation of radiometric measurements for temperature 
profiling requires the ability to accurately calculate the 
average transmittance across the sounding channel. Therefore, 
an accurate calculation of atmospheric absorption due to broad 
banded or continuum absorption features is as significant as 
one of the features due to individual lines. The important 
broad banded features affecting atmospheric absorption in the 
4. Sum region are due to the pressure induced N 2 fundamental 
{called the N 2 continuum) and to wings of relatively nearby very 
saturated CO 2 lines. 

An approximate analytical ' form for the transmittance 
through an atmospheric path, at zenith angle 0 from pressure P 
to P' (in atmospheres), due to the N 2 continuum is given by 
(Susskind and Searl, 1977) 




P"dP"] 


( 1 ) 


where T(P) is the temperature at pressure P, is the atmosphere 
scale height at 273K=7.99 km, is the mixing. ratio of N2=.79, 
and k^ is the self-induced absorption coefficient of N 2 ( amag 
km"^); The values of k.^ used, in this work are derived from 
Figure 3 of Shapiro and Gush (1966). Sample values of k^ are 
.028 at 2180 cm“^ and ,126 at 2420 cm”^. 


It is well known that a Lorentz line shape seriously 
overestimates the absorption of CO 2 lines at frequencies far 
from the line centers. Winters, Silverman, and Benedict (1964) 
and Burch (1969) have introduced empirical sub-Lorentz line shapes 
to provide an adequate fit for laboratory absorption spectra of 
CO 2 in the region near 2400 cm“^. Both line shapes can be 
expressed as 


k(v) =kj^{ v) x( I ) (2) 

where kj^(v) is the Lorentz line shape and X(Iv-VqI) is a 
multiplicative correction factor depending only on the distance 
in frequency from the line center. 

Burch presented X in graphical form, shown in Figure 1. 
Burch does not state the value of Iv-VqI above which the line 
shape becomes sub-Lorentz. The calculations in this work using the 
Burch X assume X=1 for 1 v-VqK 0.5 cm“^ . Winters, et al . 
used the following analytical form for X: 
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X I) 





( 3 ) 


^WSB ^ I ^ I ^ I K5 . 0 cm~^ 

^WSb( Iv-Vo 1)-5.) iv-v^!>5. cm-1 


The WSB X function is always larger than the Burch X and is 
also a smooth curve while the Burch X contains a relatively 
constant segment between 23 and 50 cm“^. It is shown that the 
use of a composite X function 


Xc( \^-^o^)=l 


Iv-VqK. 5 cm' 


“1 


Xc( Iv-Vol>-5 cm-1 


(4) 


having the analytical form of X^gg but agreeing more closely 
with the values of Xg/ gives substantially better agreement of 
calculated and observed spectra. 

The observed atmospheric absorption spectra shown in 
Figures 2 and 3 were taken by Beer using. an interferometer with 
0.059 cm ^ spectral resolution. The surface pressure was 793 mb 
and the solar zenith angle was 59°. More details are given in 
Susskind and Searl (1977). The theoretical calculations of 
atmospheric transmittance shown in Figures 2 and 3 were done 
according to the method of Susskind and Searl (1978),. The 
monochromatic transmittances' were convoluted with a 0.059 cm”"^ 
triangular response function. 

The main atmospheric absorption from 2385-2440 cm“^ 
arises from broad band features du.e to wings of the strong CO 2 
lines in the vicinity 2325-2375 cm”'^ and to the N 2 continuum. 
Superimposed on these are a series of weak CO 2 lines comprising 
the band head of the band. Figures 2a and 2b show. a com- 
parison of the calculated and observed spectra from 2385-2405 cm“^ 
using the Burch and WSB line shapes'. The WSB line shape produces 
systematically too much absorption over the whole region. The 
Burch line shape gives much better overall agreement but has too 
little absorption in the region 2390-2400 cm”^. This behavior 
is a result of the flat portion on the Burch X function with 
too little absorption below 23 cm“^ and too much above 50 cm“^ . 
Figure 2c, showing the results using the composite x function 
(Equation 3), gives substantially better agreement with experiment, 
though it produces slightly too much absorption below 2392 cm“'^. 

The spectra give very little indication of the value of 
beyond which the line shape becomes sub-Lorentz. A function similar 
to Xc becoming sub-Lorentz after 5 cm~^ produced only slightly 

worse agreement with experiment. 

Atmospheric absorption in the region 2180 cm“^-2240 cm“^ , 
shown in Figure 3, is due almost exclusively to lines of N 2 O, CO 2 , 
and H 2 O. The background absorption is much less than in- the 
2400 cm“'^ region and arises primarily from the N 2 continuum. The 
use of any of the sub-Lorentz line shapes in Figure 1 predicts 
negligible CO 2 absorption in this region. Agreement of theory 
with experiment is seen to be quite good. Use of a Lorentz line 
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Figure 2. Comparison of observed atmospheric absorption 
spectra (0,59 cm“l resolution) with calculated 
spectra using line shapes of Burch et al . , WSB, 
and that derived in the present work. The 
surface pressure is 793 mb. The zenith angle 
is 59". 
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15% Reduction in Width of IM2O 
















shape for CO 2 / however, or neglect of the N 2 continuum, 
results in significant errors in transmitttance calculations, 
as shown in Fig.ure 4.. 

The line parameters taken from the AFCRL compilation 
appear to be in general quite good. With the .exception of the 
H 9 O line at 2184.4 cm~^ , the calculated features agree well with 
observations. Most weak lines in the observed spectra not found 
in the calculated spectra can be attributed to CO absorption on 
the sun. The line widths of N 2 O may be too large and were reduced 
by 15 percent to get the agreement shown in Figure 3. 

While agreement between observed and calculated spectra 
is generally satisfactory, more work must be done on obtaining 
the best CO 2 line shape. For example, while Figures 1 ‘ and 2 
show a composite X which gives good agreement with observed 
spectra from 2390-2400 cm Table 1 shows two other composite 
line shapes, X 2 X 3 r which give essentially equal agreement 

of theoretical and observed spectra in the same spectral region. 
These two line shapes are almost identical to X beyond 
20 cm~'^ from the line center, but are considerably different 
less than io cm“^ from the line center. The strong CO2 lines 
whose wings give significant contribution to the absorption 
in the spectral region begin at about 2375 cm”-^. Therefore, 
calculated absorption spectra at about 2395 cm~^ using all 
three line shapes are about the same. Significant differences 
occur at lower frequencies however which are essentially 
completely opaque in the observed ' spectra; hence these spectra 
cannot be used to select an appropriate line shape. Table 1 
gives differences in clear column radiances calculated for a 
standard atmosphere using through X 4 as compared to 

the Burch X. Both the Burch x X 4 (WSB X) have been 

shown to be incorrect by Figure 2. Therefore, the relative 
error of 42 percent in calculated clear column radiance between 
X^ and Xg only shows the need .for care in performing the 
calculations. More significant however are differences between 
results using X^ through X^ which, at this time, are all 
equally reasonable line shapes. Agreement in calculated clear 
column radiances is of the order of 1 percent for the channels 
with centers 2394.5 and 2392.35. Differences become quite 
significant in the more opaque portions of the band, reaching 
9 percent at 2386.10 cm“-^., which is well above the acceptable 
level of 1 percent. Consequently, experiments are currently 
underway in analysis of high resolution balloon-borne' atmospheric 
absorption measurements which will enable us to examine spectra 
in the region 2380-2390 cm“^ for a more precise determination • 
of an accurate X function. 
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Table 1. U. S. Standard Atmosphere Percent Difference 
in Calculated Radiance Relative to Burch X. 


v{ cm“^) 

dt/dlnP MAX (mb) 

>^1 

>^2 

X 3 

X4 

2383.75 

270 

- .75 

- 3.85 

- 7.33 

-13.04 

2386.10 

390 

- 4.15 

-8.13 

-13.40 

-29.50 

2388.20 

550 

- 8.30 

-11.19 

-15.85 

-40.89 

2390.20 

700 

-10.16 

-11.58 

-14.36 

-42.48 

2392.35 

975 

- 9.50 

- 9.91 

-10.89 

-36.27 

2394.50 

SURFACE 

- 7.77 

- 7.85 

- 7.97 

-28.49 


419 

X2(X)=e“* 647(X-2) 

.390 

X3(X)=e-*749(X-5) 


• U U 

X 4 (X}-e"^^°(^-5) (WSB) 
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SIMULATION STUDIES OF THE ADVANCED METEOROLOGICAL 
TEMPERATURE SOUNDER (AMTS) 

(A, Rosenberg and J. Susskind) 

The AMTS is a new generation, high-resolution passive 
temperature sounder designed jointly by GSFC and JPL. The 
objectives of this study are to develop techniques for the 
analysis of AMTS data and investigate its performance charac- 
teristics relative to HIRS. 

The theoretical considerations and the full details of 
the AMTS are given in Kaplan, et al . (1977). It contains twelve 
temperature sounding channels in the 4 . 3vm and ISym regions, 
capable of sounding temperatures in the stratosphere and .in the 
troposphere. The spectral resolution and the weighting functions 
are shown in Figure la. It contains also three lower tropospheric 
sounding channels in the I5|im region for cloud filtering, as well, 
as water vapor sounding channels and three 4.0um window channels 
fo.r determining surface temperature, surface emissivity, and re- 
flected solar radiation (not shown in Figure 1). . For comparison, 
the weighting functions for the HIRS channels are shown in 
Figure lb. The AMTS represents an improvement over previous 
sounders such as the HIRS because of its narrow weighting functions, 
obtained by utilizing small spectral band passes in the wings of 
CO 2 lines, with high transmittance dependence on temperature and 
pressure. Consequently, it is possible to increase the number of 
sounding channels, in particular in the stratosphere, and obtain 
a higher vertical resolution. This property can be best demon- 
strated by simulated observations and retrievals based on reported 
significant level radiosonde data containing the full fine structure 
of the temperature profile. 

The general retrieval procedure used is based on the 
relaxation method outlined by Chahine (1970). In the presence of a 
cloud field, it is necessary to have observed radiances in two 
fields of view, R^ and R 2 for all the channels i. These 
"observations" are ^ simulated by calculating the radiances R^ ^ 
and R 2 for all the channels i using radiosonde temperature' 
profil4s and single layer clouds with latitude dependent 

cloud top heights and a random distribution of fractional cloud 
covers, N^ and N 2 , in the two fields of view. A Gaussian distribu- 
tion of random noise, representative of the performance of the HIRS 
and AMTS instruments, is added to the calculated radiances of the 
respective channels. 

From the observed radiances in the two fields of view, 
clear column radiance Rj^ can be reconstructed from the relation 


^i " %,i ^ ^^^l,i ■ ^2,i^ 


( 1 ) 
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Figure 1. Normalized temperature weighting functions for HIRS and AMTS for nadir viewing. 






The problem of temperature retrievals in cloudy 
atmospheres involves the simultaneous determination of n from 
Equation (2) and a temperature profile, {Tj^}, such that 
Equation (1) will prevail within noise levels for the same 
value of n for all the channels. h are obtained by 

using the following iterative method. For the iteration, the 

clear column radiances, are calculated from the n tem- 

perature profile {T } for the temperature sounding channels (t) 
and for the cloud filtering channels (c). is then determined 

from the cloud filtering channels by using in Equation (2). 

Using Equation {1) and again, the reconstructed 

clear column radiance is determined for each temperature sounding 
channel. These values are then compared with the clear column 
radiances computed from the n^“ temperature profile. If 


e 


l^t 


(n) 



(n) 


<.001 , 


the n^^ temperature profile is taken as the retrieved temperature 
profile {Tj^} = {Tj^ }. If £>.001 the iteration is continued 
by constructing a new temperature profile, from {T } 

in the following manner. The temper atur es T(P^t corresponding to the 
pressures, of the peak of the weighting functions for channel t 
are corrected according to the equation : 




(3) 


The full 66-level temperature profile used to calculate , 

{Tj^^q^}, is obtained by interpolating between the points T(P^) 
linear in the logarithm of the pressure, while preserving the shape 
of the first guess temperature profile (see next section). If 
£>.001 after ten iterations, the 10th temperature profile is taken 
as the retrieved profile {Tj^} = {T^g}. 

The iteration process must start with a first guess 
temperature profile (FGTP) {T^}. Since the shape of {T^} is 
pre-served during the iteration procedure, it is important that 
the shape of {T^} will be as close as possible to that of 
the true profile for a good solution to be obtained. 

.This requirement is particularly important in the tropopause 
region which is characterized by a sharp discontinuity in the 
lapse rate -{3T/9Z). 
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For HIRS, a zonally-averaged climatology profile was 
used as an FGTP. The deficiency of this FGTP is that it does not 
reflect the shape of the individual radiosonde temperature profiles. 
For the AMTS, an innovative method was developed, making use of the 
sharpness of the weighting functions to construct an initial FGTP 
that does not require- any prior knowledge of the state of the 
atmosphere . 


The key to the method is the assignment of a characteristic 
pressure, to each channel, representative of a pressure whose 

temperature has been shown to be close to the observed brightness 
temperature by examination of results from a variety of profiles. 

These characteristic pressures, do not coincide with P^, 

the pressures at the peak of the weighting functions. The first 
guess is constructed by assigning the observed brightness tempera- 
tures, T^, to the corresponding characteristic pressures 
The full 66-level temperature profile is then obtained by a linear 
interpolation in the logarithm of the pressure. 

This temperature profile is further refined to reflect 
the shape of the true profile in the tropopause region. Starting 
from the surface, the lapse rate is calculated for each segment, 

P^^^^ to until a segment L is reached in which the lapse 

rate is less than 2®K/km, The temperature profile in this segment 
is modified by replacing it with two linear segments, which are 
extrapolations of the two adjacent segment L-1 and L+1 (see Figure 2). 
The values of the temperature and the pressure at the intersection 
of the segments (L-1) and (L+1) are in general fairly close to the 
corresponding values of the true tropopause. 

The FGTP obtained in this manner for the AMTS channels 
together with the triangulation is greatly superior to the latitude- 
dependent climatology profile. 

This method of constructing an FGTP and obtaining approxi- 
mately the height of the tropopause by triangulation was not applied 
to the HIRS because of the broadness of the weighting functions of 
its stratospheric sounding channels. Therefore, a climatology 
temperature profile was used as an FGTP. However, studies are 
underway co modify this method to test whether it is applicable to 
the HIRS as well. 

With the exception of the choice for the FGTP, the retrieval 
procedure for the two sounders was identical. After a retrieved 
temperature profile {Tj^} was obtained, the RMS errors relative 
to the profile {Tr^} were calculated for the 12 mandatory levels. 

The overall and the mandatory levels of RMS errors obtained from 286 
global radiosonde stations are given for the AMTS and for. the HIRS in 
Figure 3. The overall RMS errors for AMTS and HIRS are 1.72°K and 
2.27®K, respectively. The AMTS is superior to the HIRS particularly 
in the stratosphere. A typical retrieved temperature profile for 
the AMTS and the HIRS together with the radiosonde profile used to 
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Construction of the first guess temperature profile using 
characteristic pressures, observed brightness temperatures, 
linear interpolation of temperature in ^nP and triangulation 
in the region of a break in the apparent lapse rate. The 
first guess temperature pressure is shown as P, . 
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lUys errors versus height of simulated AMts 
and HIRS retrievals as compared to radiosonde. 





simulate the observations are given in Figures 4a and 4b, 
respectively. For the AMTS the true and the retrieved tropo- 
pause are also given. That, portion of the AMTS profile in 
which triangulation was used to determine thre tropopause is 
shown by do'ts. As is expected, it can be seen that the AMTS 
more accurately reproduces the structure of the radiosonde 
profile, particularly in the tropopause region and the lower 
troposphere, due to the sharpness of its sounding channels. 

The RMS error for the retrieved tropopause heights 
was calculated for the AMTS and was found to be less than 
1.0 km. The general accuracy of this determination can be seen 
in Figures 5a and 5b showing contour maps for the tropopause 
pressure over the U. S. as determined from radiosonde reports 
and simulated AMTS retrievals. The general agreement is seen 
to be very good . 
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Figure 4. Sample radiosonde profile compared to retrieved AMTS and KIRS profiles. The 
dotted portion of the retrieved AMTS profile v:as determined by triangulation. 










RAPID ALGORITHM FOR CALCULATION OF TEMPERATURE 
DEPENDENCE OF HIRS TRANSMITTANCE FUNCTIONS 

(T. Mo and J. Susskind) 

The purpose of this study is to establish a rapid 
and accurate algorithm for calculating the temperature and 
zenith angle dependence of the dry atmospheric transmittajice 
functions for 11 HIRS channels in the Nimbus-6 satellite. 

Analytic formulas of atmospheric transmittance func- 
tions for 11 HIRS channels are developed and parameterized 
as functions of atmospheric temperature profiles. Parameters 
for atmospheric transmittance functions "at two solaf“zenith 
angles of 0° and 40° are obtained from least-square fit to 
line-by-line calculated transmittance functions. Based on the 
two sets of parameters at the angles of 0° and 40°, transmittance 
at an arbitrary zenith angle 6 is obtained by linear inter- 
polation in sec0. Atmospheric radiances calculated with the 
transmittances obtained from the analytic formulas show small 
errors of the order of 0.5 percent compared to line-by-line 
calculations. These analytic formulas have been incorporated 
into the analysis of HIRS observations for determination of 
atmospheric temperature profiles. 

The dry transmittance function XiCP^f 0 ^ T, 6) , from 
the top of the atmosphere to pressure at zenith angle 6, 
averaged over channel i, is a function of the entire temperature 
profile from pressure Pj^ to the top of the atmosphere. For 
modeling the temperature dependence of the transmittance function, 
it is convenient to introduce an effective layer transmittance 
for a layer between pressure P^ and P^-i defined by. 


TjL(Pn/Pn-l^'^' = 'Ti{Pn^0/T, 0)/T^(Pn_j^,O,T, 0) 


( 1 ) 


For simplicity, this effective layer -transmittance is written as 

Ti,n(T,e). 

For monochromatic radiation, the effective atmospheric 
transmittance for the layer is identical to the atmospheric 
transmittance through the layer. It therefore depends only on' 
the temperature profile in the layer (taken to be homogeneous 
for thin layers) and is independent of the temperatures above 
the layer. While this property is not rigorously true for 
finite width channels as in the HIRS case, the effective layer 
transmittance is still dependent primarily only on the tempera- 
ture of that layer with only a small residual dependence of 
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the temperature profile above the layer. The procedure used 
to model the temperature dependence of the HIRS transmittance 
functions was to calculate transmittance func.tions- for a number 
of climatological proflTes using an accurate line-by-line trans- 
mittance program, and model the calculated layer transmittances 
as a function of layer temperature and another quantity, which 
may be considered as a mean lapse rate above the layer. A 
66-layer pressure mesh was used- in all calculations and the 
effective transmittances were calculated using adjacent layers. 

A.. Line-by-Line Calculation. The calculations of 
T^(p ,O,T,0) are done according to the method of Susskind and 
Searl (1978) using the 1976 version of the AFCRL compilation of 
atmospheric molecular line parameters. For each channel, the 
transmittance function was calculated with six climatological 
temperature profiles corresponding to the U. S. standard atmos- 
phere, January 40®N, 70'*N and July 0°N, 40°N, 70°N. Each 
calculation was done at two zenith- angles 0° and 40 “N, 
corresponding to the extreme angles of the ,HIRS instrument. 

A Voigt absorption line shape is assumed for all 
atmospheric species except CO 2 . Following a recent study by 
Susskind and Mo (1978), an 'empirical' CO 2 line shape is adopted 
and defined by 


k(v) = k^(v) X(A) 


( 2 ) 


where kj^(v) is the Lorentz line shape and X ( A)=X ( | v-v^ | ) is an 
empirical correction factor defined_as 


X(A) 


^ ^ *1 
1 for A^O.5 cm ^ 

^-a(A-0.5)^ for A>0.5 cm 


where a = 0.601 and 8 = 0.434 are two empirical parameters, 
determined to give the best fit to an observed atmospheric 
absorption spectra in the wavenumber region, of 2388-2395 cm“ 
(see Susskind and Mo, 1978). The Voigt function is used for 
A<.5 cm 


(3) 


B. Analytic Representation . The analytic formula 
employed in the present work to represent the line-by-line calcu- 
lated effective layer transmittance functions, T.: (T, 9) , has 

the form ■ ' 
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(4) 


Ti^n(T,0) = 

where the three parameters cr^ “i 

independent of temperature anc5 are obtained by least-square fit 
to the line-by-line calculated layer transmittance functions for 
each channel at each of the two angles. The quantities Aj^T and 


A^T are defined as ' 


A T 
n 




(5) 


and 



A T 
n 


A T 
n-m 


( 6 ) 


where T is the temperature at the mid-level of the layer n, and 
Tg corresponds' to the U. S. standard temperature at the layer n. 

The quantity in Equation (6) is taken to represent the difference in 
mean lapse rate of the profile at pressure Pj^ as compared to that 
of the standard atmosphere. This quantity was difficult to model 
accurately but the following choice of m values gave reasonable 
fits for a variety of profiles. 


m = 


0 

n-25 


i f n;<l 0 , 
if 10 <n_<3 0 , 
if n>30. 


(P<10 mb) 

(10 mb<P^200 mb) (7) 

(P>200 mb) 


The first two terms in the expression on the right hand of 
Equation (4) are essentially a Taylor series expansion of the 
layer transmittance at a pressure level P, using the temperature 
T as variable. Approaches similar to this two- term approximation 
have been employed by Kondrat'ev and Timofeyev (1969), and by 
Wang (1978). However, the last term in Equation (4) was found to 
be important in obtaining substantial imprpvement in the agreement 
between the line-by-line calculations and the analytic results. 

Physically, the last term in Equation (4) takes into 
account the temperature 'memory' -effect of the upper atmosphere 
on the transmittances calculated at lower layers. Its importance 
becomes more pronounced’ as the pressure increases in the lower 
atmosphere. Equation (7) shows that this term (with m=0) has no 
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effect on the transmittances in the upper atmosphere where P^lO'mb 
(i.e., the effective lapse rate is always taken as zero in the 
upp^r stratosphere) . The effective lapse rate correction is taken 
over a 5-50 mb interval for pressures between 15 mb and- 200 mb, 
and taken over the interval between P and 150 mb for P>200 mb. 

Equation (4) was used to fit the line-by-line calculated 
transmittance at each of the 66 levels (i.e., n=l-66), and a set 
of best-fit parameters for cf^^, and was obtained for each 
level n, channel.! and angle 0. The layered transmittance of an 
arbitrary angle 9 between 0® and 40® is obtained by linear inter- 
polation of the values of In t in secS. 


Given the layered transmittance at any angle, the 
total atmospheric' transmittance up to a pressure level P can 
be obtained from Equation (1), 


T^(P^,0,T,e) = (8) 


The agreements between the results obtained from 
Equations (4) and "(8) using the best-fit parameters, and the line- 
calculated total transmittance as function of pressure is shown 
in Figure 1, In this figure, the solid curve represents the 
total transmittance calculated with- the best-fit parameters, and 
the dots are obtained by line-by-line calculations at 6 = 20°. 

The total errors (including the interpolation in sec9 ) are of 
the order of 0.01 percent and are much smaller than the sizes of 
the dots. The temperature profile used in this calculation is 
taken from Rodgers and Walshaw (r9'65) and corresponds to. the 
actual temperature of Nairobi, January 7, 1960, 


To demonstrate the accuracy of the fitted results and 
the usefulness of the analytic formulas developed in this study, 
the total atmospheric radiances for the 11 HIRS channels^ V7ere 
calculatei.; by using the total transmittance from. the exact line-byr- 
line calculation and analytic formulas (see Equation 8). The 
comparison of these calculations is shown in Table 1, which a.lso 
gives the percentage errors of the calculated radiance resulting 
•from the analytic formula. It can be seen from this table that 
the magnitude of the errors is less than 1 percent. The last 
column in Table 1 gives the largest errors, which occur among the 
66 total transmittances t (P j^, O,T,0 ) for each channel i if the 
analytic formula is employed. In practical application, the 
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Figure 1. Comparison of the line-by-line calculated 
transmittance with 'that obtained from 
best-fit parameters. 
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Table 1. Comparison of the total atmospheric radiances obtained 
from the exact line-by-line calculated transmittance.,, and from 
the analytic formulas developed i-n this study. The temperature 
.profile used in these calculations is taken from the data given 
by Rodgers and Walshaw (1966), and corresponds to the temperature 
profile of Nairobi, January 7, 1960. The zenith angle of observation 
is 9 = 20°. The error (in percent) includes that of interpolation 
in angle for obtaining transmittances at 9 = 20°. 


^0 

( cm”^ ) 

Exact 

Radiance 

(erg sec“^cm”^s-r“-’-) 

Calculated 

Radiance 

(erg sec“^cra“^sr“^ ) 

Error 

(%) 

Largest 

Error 

in 

669.0 

58.173 

58.153 

0.03 

-. 002 

679.0 

40.328 

40.341 

-0.03 

-.001 

690.0 

40.413 

40.461 

-0.12 

-.002 

700.0 

50.876 

51.123 

-0.49 

-.005 

716.5 

73.626 

74 ..074 

-0.61 

-.006 

732.0 

90.287 

90.612 

-0.36 

-.003 

749.5 

109.654 

110.017 

-0.33 

-.004 

2190.0 

2-. 009 

2.014 

-0.25 

-.002 

2210.5 

1.138 

1.141 

-0.22 

-.001 

2243.5 

0.400 

0.402 

-0-.69 

-. 002 

2357.0 

0.128 

0.128 

-0.04 

-.001 
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instrumental measurement of radiance usually has an accuracy of 
about 1 percent. The analytic formula gives results with an 
accuracy well' within the instrument-error limit, and has been 
incorporated in the analysis of HIRS data. Use of this new 
analytical form has produced a significant improvement in the 
accuracy of retrieved temperature profiles and has removed the need 
for empirical tuning of transmittance functions. 
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DETERMINATION OP ATMOSPHERIC WATER CONTENT AND 
ITS EFFECT ON SCAMS TEMPERATURE RETRIEVALS 

(J. Rosenfield and J. Susskind) 

The present, temperature retrieval methods employed at the 
Modeling and Simulation Facility make use of the observations from 
all five microwave channels of the Nimbus -6 SCAMS instrument. 
Channels 1 and 2 , at 22.235 GHz and at 31.4 GHz, are window chan- 
nels whose observations are primarily determined by the surface 
temperature, surface emi-ssivity, and atmospheric liquid water' and 
water vapor content. Channels 3, 4, and 5, at 52.85, 53.85, and 
55.45 GHz, lie in the wing' of the 60 GHz band of oxygen and sound 
atmospheric emission from roughly 950 mb, 500 mb, and 200 mb, 
respectively. While brightness temperatures of these channels are 
primarily a function of the atmospheric temperature profile, the 
observations in channel 3, and to a lesser extent channel 4, are 
also affected by surface emissivity and atmospheric water content. 

In the current operational procedure, these quantities are deter- 
mined from the brightness temperatures of channels 1 and 2 , 
using an empirical scheme of Waters, et al . (1975). This scheme 
ignores the effects of atmospheric liquid water content on the 
observations. This study was undertaken to extend the method of 
Waters, et al . , to' include the effects of atmospheric liquid 
water content on the SCAMS observations. In addition, a technique 
was developed for the determination of atmospheric liquid water 
content and water vapor content from these observations. 

■Theoretical calculations of water vapor transmission 
in the microwave region were performed for 124 temperature-humidity 
profiles. The results verify the following empirical relationship 
of Waters, et al ., between the water vapor plus oxygen transmissions 
in the Nimbus -6 SCAMS microwave channels 1 and 2, 

T 2 (v, 02) = .32 Tj (v,02) + . 666 , (1) - 

to within +.005 for nadir viewing, with slight degradation at 
zenith angles up to 50 The following linear relationship 
between the water vapor only transmissions in channels 1 and 3 
was found _o be accurate to .007: 

T 3 (v) = .55 (v) + .45 (2) 

In addition, the channel 1-3 water vapor transmissions were modeled 
as a function of the water vapor column density, W, as follows: 


-a.W 

T.(v)=e ;i = 1,2,3 
at = .065 
a, = .019 
33 = .036 

where W is in g/cm^. 


(3) 
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In orde^r to evaluate the effect of clouds on the current 
temperature retrieval scheme, theoretical calculations of brightness 
temperatures, Tg, for channels' 1-3, with and without the inclusion 
of clouds, V7ere performed. The absorption arising from cloud water 
droplets was included, but scattering effects, which become signi- 
ficant only in the case of heavy rain, were neglected. A single 
cloud layer was placed between 1 and 2 km (~ 850 mb) with water 
droplet densities, pg, of .1, .2, .3, and .4 g/m^ , corresponding 
to liquid water contents, L, of 10, 20, 30, and 40 mg/cm^ . Over 
ocean, brightness temperatures in channel 3 increased by ~2-3°K 
with clouds corresponding to L up to 40 mg/cm^. Over land, clouds 
affect the channel 3 brightness temperature by only .1-.2°K because 
land emissivities are close to unity. 


The extrapolation techniques of the current temperature 
retrieval program were tested to determine the errors incurred 
when clouds are not taken into account. Over ocean, when no 
clouds are present, reconstructed channel 3 brightness temperatures, 
Tb 3 's, are -.1• ** different from exactly calculated Tg^'s. Vlith 
clouds, errors in Tg^ sre of the order of 1°K for L = 30-40 mg/cra^ . 
Over land, errors in Tg^ are ~.1-.2®K, with or without clouds. 

In order to modify the current operational temperature 
retrieval program to take into account the effects of dense clouds 
on the SCAMS channel 3 brightness temperatures, a procedure has 
been developed to determine estimates of water vapor and liquid 
water content from observed brightness temperatures of the SCAMS 
microwave channels 1 and 2. As a prerequisite to this, it was 
necessary to model the. following three quantities: 


• 1. Fith a single cloud layer 
emission temperatures, T^^ and T^^ 2 ^ 
mately given by ' ' 


, the effective atmospheric 
were found, to be approx i- 


'^AT,1 ~ ^s' ^AT,2 “ -s' (4) 

where Tg is the surface temperature. 

2, The following empirical relationship between the 
liquid water transmissions in channels 1 and 2 was found: 

T2(^-) = 1.92 {H) - 0.92 (5) ' 

3. Liquid water transmissions could be modeled by 


exp[-b^L exp{ 2140/T^) /T^] 
v/here b^^ = .00012, b 2 = .00023, 
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= cloud tenperature, and L = liquid water content in ng/cn^ 


Given the neasurenents Tg^ and '^ 02 / 
and T 2 (v, 02 ,. 2 ') can be reconstructed as follows: 

Tia,v,G2) = ^V^AT,i^ ^ 




{(' 


AT 


^i-T3)2(l-Ri)2 _ 4T^T,i^i('^Pi-'^AT,i)"> /2TAT,i^^i^ (7) 


where R* is the reflectivity in channel i which, over ocean, can 
be modeled accurately given an estimate of the wind speed. Then, 
knowing that 

Ti(v,02,5,) = (8) 

and from (1) and (5) one can solve for Tj|^(v, 02) in terms of 
the reconstructed t^(v, 02 ,^-) and T 2 (v, 02 /i) as follows: 

'Tj(v, 02) = 1.0435 X3^(v,02,Jt) - 1.6984 T2 (v,02,A) - 1.0406 

+ 1.6984 {[.61^4 Tj^(v,02^JI) - T2(v,02,JU - .61272]2 

+ 1.50582 Tj(v,02,A)^} ■ 

Prom Equation (8) it follows that 



Tj^(£) = Ti(V,02,Jl)/Ti(V,02) . 

(10) 

also 

Ti(v) = t^(v,02)/t^{02) , 

< 

where 

Xj^( 02 ) has been modeled. 


water 

using 

From the quantities x^(v) and x^(£) 
content, L, and water vapor content, W, 
Equations (3) and (6). 

the liquid 

can be reconstructed 


The agreement of the reconstructed water vapor and 
liquid water content with the model is very good over ocean. 

RMS errors 'in L are 1.61 mg/cm^ , and in W are 0.045 g/cm^ . This 
agreement is an order of magnitude better than that obtained 
with schemes (Staelin, et al . , 1976; Grody, 1976) using 
statistical methods. Over land, the error in W is "1 g/cm^ 
and in L is ~ 12 mg/cm"^. These errors are comparable to 
those obtained by statistical methods. This is a result of the 
fact that the effects of liquid water on SCAMS observations over 
land are small. This being the case, an accurate determination 
of the liquid water content ovex land is unnecessary for temper.aT-. 
ture sounding. 
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The sensitivity of the reconstructed W and L to errors 
in the estimated values of the surface temperature, Tg, the cloud 
temperature, T^, and the emissivity, e, were examined. An error 
in Tg of 3°K leads to ~5 percent errors in W, ~10 percent 
errors in L; an error in of S^K leads to ~15 percent errors 
in L; and an error in e of .01 leads to ~10 percent errors in 
W, "'10-30 percent errors in L. 
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IMPROVED CALCULATION OP SCAMS BRIGHTNESS TEMPERATURES 


(S. Kirschner and J. Susskind) 

The SCAMS tenperature-humidity sounder on Nimbus-6 is 
a five-channel passive microwave radiometer. For optimized 
retrieval of atmospheric temperature profiles from SCAMS 
measurements, it is necessary to be able to accurately calculate 
the brightness temperatures for a channel given the temperature 
profile. The tirst two channels on SCAMS are atmospheric 
windows and their radiances are influenced primarily by the 
surface temperature, surface emissivity, and water vapor and 
liquid water content of the atmosphere. The remaining three 
channels are in the O 2 absorption band and sound thermal emission 
from the earth's atmosphere, beijig_sensitive jp^imarily to tempera- 
tures near the peaks of their weighting functions. The' channel 
frequencies are shown in Table 1 together with the peaks of the 
weighting functions and surface transmittances- for nadir viewing. 

Channel 3, and to a lesser extent Channel 4, see the 
surface, and these brightness temperatures are therefore affected 
by surface emissivity as well as atmospheric temperature profile. 
In our analysis of SCAMS data, the surface emissivity is estimated 
from the observations of Channels 1 and 2, and a guess surface 
temperature, using a procedure developed by Waters, et al . (1975). 
This technique assuming atmospheric attenuation of radiation in 
Channels 1 and 2 is due only to water vapor and O 2 , and ignores 
the effects of liquid water clouds. The two observations are 
used to determine the surface emissivity and atmospheric trans- 
mittance in Channel 1, and the corresponding quantities in 
Channels 2 and 3 are related to those in Channel 1 by theoretical 
expressions. 

In order to obtain accurate temperature profiles from the 
SCAMS observations, one must be able to calculate brightness 
temperatures as a function of temperature to an accuracy of 
about ,5°C. This involves a knowledge not only of the surface 
emissivity but also accurate atmospheric transmission functions 
for the channels. The attenuation due to O 2 absorption is calcu- 
lated accurately according to the theory of "^Rosenkrantz (1975). 
Additional atmospheric attenuation arises from absorption by 
liquid water and water vapor. Both effects are small and it was 
originally thought that the common approximation of neglecting 
these effects in the calculation of brightness temperatures for 
the temperature sounding channels was valid, i.e., that the 
differences in calculated radiances for a given temperature 
profile, including and excluding the effects of atmospheric 
water, was less than ,5®C for Channels 3-5. In accordance with 
this, attenuation due to water vapor was ignored in the temper- 
ature retrievals done for the DST 5 and 6 periods. 

Comparison of retrieved temperatures with colocated 
radiosondes (mostly over land) produced good results and demon- 
strated that the approximation was satisfactory. A further study 
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Table 1. ■ SCAflS Sounding Channels. 


Channel 

Frequency 

(GHZ) 

Peak of 
Wt. Function 
(mb) 

Surface 

Transmittance 

1 

22.235 

window 


2 

31.400 

window 


3 

52.850 

950 

.33 

4 

53.850 

500 

.10 

5 

55.450 

200 

.002 
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of retrieved surface temperature over the tropical Pacific showed 
large systematic differences between temperatures . retrieved using 
SCAMS and HIRS (infrared) observations. This l.ed to the examina- 
tion of the effect of the approximation more closely, and it was 
found to be extremely bad over tropical oceans but good over land, 
as shown in the next section. 

Sensitivity of Brightness Temperatures to Water Vapor 

The brightness temperature observed by a downward looking 
microwave passive radiometer sounding a frequency v- can be 
written as - - i 


Ti = £iT3Xi (P^ }+ ( 1- e . ) (1- Ti (P^) )Fi ^T^ 


(Pg)+(1-T. (Pg) )F.Tg 


( 1 ) 


where the first term represents the radiance being emitted by 
the surface, with emissivity temperature T , and atmospheric 
transmittance Tj^(Pg), the second term represents downward 
atmospheric flux reflected from the suface, and the last term 
represents radiation emitted by the atmosphere. The terms 
^i^s ^i^'^s equivalent atmospheric temperatures 

looking down from the satellite and up from surface, respectively, 
and represent a weighted mean temperature over the region of 
atmospheric emission. In the limit of an optically -thin atmosphere 
such as found in Channels 1 and 2, F^ = F^-*. Otherwise, F- ■- is 
greater than F-, that is, the mean emission of radiation in the 
downward direction takes place lower in the atmosphere than that 
in the upward direction. 

The atmospheric transmittance in Equation (1 ),t 4(P ), 
differs from "1" due to attenuation by O 2 and H^O vapor (liquid 
H 2 O drops are being neglected in this discussion) and can be 
written as a product 


= hi,ox(P) (2) 

where Tj^^Qjj(P) and represent the nadir transmittances 

due to oxygen and water vapor, and 0 is the zenith angle of 
observation, x calculated as a function of temperature 

profile according to Rosencraritz (1975). Given x* (P )/ x- (P) 
is calculated according to the approximation of a^SVl ^ wa^4r 
vapor scale height and temperature independent absorption 
coefficients/ i.e. 


i ,w^^^ ~ i ,w 


(Ps) 


4.8 


(3) 


It is clear that both addition of water vapor and increase of 
zenith angle of observation decrease the surface transmittance 
Ti(P_) and therefore modify the observation as shown in 
Equation Jl). A more subtle effect is the small dependence of 
Fj^ and F^'' on water vapor and zenith angle. Increasing 0 


/*• 
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will increase the optical path of the atmosphere homogeneously 
and therefore increase F-^ and decrease P^. Increasing the 
water vapor amount also increases the optical path of the 
atmosphere, but primarily at lower altitudes. Therefore, both 
and increase with increasing water vapor. 


A change in brightness temperature with t (addition 
of water vapor) can be written as 

dT. 

-±~— - (G ,-H . )T 

dT^(-Pg) ^ 1 i s 


where 

Gi = £i+(l-ei) 

and 


(l-2Ti(Pg)F{+Tj^(Pg) (l-Ti(Pg) 


dTi(P^) 


(4a) 


(4b) 


Hi = P. - (l-Ti(P3))^ 


dx. (P ) 
1 s 


(4c) 


The two terms in Equation (4a) , and H- represent the 
increase in radiation coming from the surface and the decrease 
in radiation coming from the atmosphere, with increasing x. 


Table 2 shows observed and calculated brightness tempera- 
tures for a portion of a scan' of the SCAMS instrument over the 
tropical Pacific. The calculations are shov/n for two models, one 
neglecting the effect of water vapor on the transmittance and 
one using water vapor determined from the observations in Chan- 
nels 1 and 2., by way of Equation (2) and the empirical equation 


^3,w 


(Pg) 


.55 



+ .45 


(5) 


The emissivity £3 was also determined from the observations in 
Channels 1 and 2. Its zenith angle dependence shows good agreement 
with theoretical values for sea water with a wind speed of about 
12 meters/sec. The temperature profiles used in the calculations 
were 12-hour forecasts from the GLAS GCM. The forecast surface 
temperature is included in the table as well as the retrieved 
surface temperature obtained from inverting the SCAMS observations. 
Also shown in the table are F, dF,. f ' , dP", G, and H for Channel 3. 

dx dx 


It is clear from the table that addition of water' vapor 
over ocean causes a large increase in the brightness temperature 
of Channel 3, which decreases with increasing zenith angle from 
2.6°K at 11 degrees to .9°K at 49 degrees. The effect is due 
primarily to the low surface emissivity over ocean and the decrease 
in surface transmittance with increasing zenith angle. The bright- 
ness temperatures calculated from the initial guess, including 
water vapor, are in much better agreement with observations than 
those calculated ignoring the water vapor absorption. While this 
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Table 2 


. Quantities calculated over water for various senith angles and water vapor 
models. The temperature retrievals are over the Pacific Ocean for 
22Z 21 AUG 75 at -2* (±2®) latitude and 152.5® (±10°) longitude west of 
Greenwich. 


Model 

Zenith 

Angle 

^3 

■^3 

T. 

3«) 

T ■ 
S 

(K) 

F 

s 

^3 

F ’ 

-dF3/dT3 

.dP3/ffT3 



(deg) 



Obs. 

Calc'.’ 

Fcst . . 

Ret’d 

(mb) 

6 

3 

3 

(i) 

11 ■ 

.516 

.321 

262.0 

259.0 

298.9 

304.1 

1013 

. 938 

.886 

.157 

.094 

.661 

.950 

(No 

24 

.515 

.294 

262.9 

260.2 

298.8 

302.6 

1015 

.941 

.885 

.164 

.093 

.686 

.950 

water 

36 

.518 

.253 

263.6 

261.6 

299.5 

301.9 

1013 

.939 

.879 

.176 

. 087 

.725 

.944 

vapor) 

49 

.544 

.181 

263.2 

262.1 

298.6 

298.2 

1016 

.946 

.873 

.209 

.083 

.805 

.941 

(ii) 

11 ■ 

.516 

.290 

262.0 

261.5 

298.9 

298.0 

1013 

.943 

.889 

.157 

.094 

.691 

.956 

(Water 

24- 

.515 

.263 

262.9 

262.5 

298.8 

297.7 

1015 

.946 

.888 

.164 

.093 

.716 

.956 

vapor ) 

36 

.518 

.222 

263. 6 

263.5 

299.5 

297.4 

1013 

.945' 

.882 

.176 

.087 

.755 

.949 


49 

.544 

.157 

263.2 

262.9 

* 



298.6 

295.3 

1016 

.951 

.875 

.209 

.083 

.828 

.945 




agreement may be fortuitous depending on the accuracy of the 
forecast guess, it is significant to note that the differences 
between observed and calculated brightness temperatures have 
negligible zenith angle dependence when water is included, but 
have a large zenith angle dependence when dry transmittances 
are used. 

Table 3 shows the same type of information for a scan 
over land, which is characterized by a- surface enissivity close 
to 1. In this case the brightness temperature calculated with 
and without water vapor differ by no more than .2°K. H, the 
dependence of atmospheric radiation on water vapor, has changed 
only slightly, but G, the dependence of the surface radiation 
coming from, water vapor, has risen substantially to be slightly 
larger than H, producing in. fact a small .drop in brightness 
temp.erature with inclusion of water vapor. This drop would 
have been substantially larger were it not for the term dF 2 ''/dT 2 i^ 
which produces an increase in the effective atmospheric 
temperature for upwelling radiation with increased atmospheric 
water vapor. 

Incorporation of this algorithm into the HIRS/SCAMS 
temperature retrieval package led to substantial improvement in 
the quality of retrieved temperature fields. 
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Table 3. Quantities calculated over land for various zenith angles and water vapor 

models. The temperature retrievals are over the United States for 08Z 19 AUG 75 
at 50° (±2®) latitude and 112.5° (±10°) longitude west of Greenwich. 


Model 

Zenith 

Angle 

(deg) 

"3 

”^3 

T 

3<K) 

T 

s 

(K) 

^s 

(mb) 

F 

3 

F ' 

^3 

-dFi/dXs 


S 

«3 

Obs. 

Calc . 

Fcst 

Ret’d 

(i) 

6 

.953 

.419 

264.4 

264.7 

282.8 

282.8 

875 

.946 

.905 


.095 

.959 

.960 

(No’ 

17 

.953 ■ 

.426 

264.6 

•264.2 

284.7 

285.1 

847 

.931 

.891 


.095 

.958 

.945 

water 

' 32 

.954 

.378 

260.6 

262.1 

284.0 

280.2 

853 

.932 

.888 

.150 

.090 

.963 

.944 

vapor ) 

45 

.956 

.277 

259.1 

260.0 

283.1 

279.5 

898 

.947 

.892 

.170 

.084 

.973 

.953 

(ii) 

6 

.953 

.381 

264.4 

264.7 

282.8 

282.5 

875 

.951 

.908 

.136 

.095 

.962 

.968 

(Water 

17 

.953 

.387 

264.6 

264.1 

284.7 

285.4 

847 

.937 

.894 

.142 

.095 

.962 

.952 

vapor) 

■ 32 

.954 

.342 

260.6 

261.9 

284.0 

280.2 

853 

.937 

.891 

.150 . 

.090 

.966 

.950 


45 

.956 

.255 

259.1 

259.9 

283.1 

279.6 

898 

.951 

.894 

.170 

.084 

.975 

.957 




























III. OCEAN/AIR INTERACTIONS 



■realistic simulations of the global observing system 

AND of SEASAT-A MARINE WIND DATA 
(M. Cane, M. Halen, V. Cardone, I. Halberstan, and J. Ulrich) 

A 3d-day history run made with the GMSF general circu- 
lation model was used to fabricate simulated observations at the 
times and locations of the conventional surface, radiosonde and 
ship reports actually received during February 1976. The 
fabricated observations, suitably degraded for instrument and 
sampling errors, were then used to create analyzed fields on the 
GCM grid in an analysis forecast cycle like those in use at 
major meteorological centers. The control fields so produced 
are much more representative of actual analyses than those pro- 
duced by perturbation of initial states with random errors. 
Significantly, the forecast error growth in five simulated 
72-hour forecasts from the control states is similar to that 
found in operational numerical forecasts. 

• Further experiments simulated the addition of surface 
winds derived from Seasat-A to the control run. The Seasat wind 
data were fabricated directly for the GCM gridpoints intercepted 
by a Seasat-A scatterome ter swath and were assumed to be error-free 
and representative of the lowest active level of the GCM (level 9). 
Two different asynoptic assimilation methods were used; direct 
insertion and successive correction. It was found that assimilation 
of the wind data by the DIM led to 10-20 percent reduction of error 
in the specification of the low-level wind field over the ocean 
relative to the control analysis, and only slight improvement in 
the simulated 3-day forecasts of sea level pressure in extratropical 
oceanic regions. Assimilation of wind data by the'SCM resulted in 
a substantial beneficial impact of the satellite data, on analyses 
and forecasts of surface pressure, and low and mid tropospheric 
wind fields in both land and oceanic regions. 
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Figure 1. Growth of the level 9 zonal wind errors (RMS) averaged over five 
simulated forecasts made from Control (C) , Seasat direct 
insertion method (DIM) , and Seasat successive correction method 
(SCM) initial states. The superior analysis (i.e., day 0) in 
the SCM case results in an improved forecast. 
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Figure 2. Growth of the sea level pressure errors (RMS) averaged over five 
simulated forecasts made from Control (C) , Seasat DIM, and 
Seasat SCM initial states. The analysis (i.e., day 0) in the 
SCM case results in an improved forecast. 










ADDITIONAL SEASAT SIMULATION EXPERIMENTS 
(M. Cane, V. Cardone', and D. Edelinann) 

This brief report updates our detailed account of 
earlier experiments simulating the global observing system and 
the Seasat-derived marine winds (Cane, et al . , 1978, submitted 
to Mon. Wea. Rev , and abstracted . in this research review). This 
work is currently being extended in two directions. 

The first is to account for the errors in the Seasat- 
derived surface winds. In a recently completed experiment, 
random errors of 2 m/s in speed and 20° in direction were added 
to the inserted level 9 winds. The winds were assimilated by 
the successive correction method (SCM). Quantitative measures 
such as RMS sea level pressure errors showed that the resulting 
forecasts were not significantly different in quality from 
those made with error-free simulated winds. Table 1 provides 
measures of 3-day forecast root mean square (RMS) errors from 
the error-free (SCM) and random error" (SCME) initial states. 

This lack of difference is attributable to the 
SCM's ability to smooth random errors. Current work is 
focused in implementing the actual scatterometer geophysical 
retrieval algorithms in order to correctly capture the unusual 
error structure to be expected from this instrument. At the 
same time the boundary layer parameterization in the GCM is 
being modified to account for the turning of the wind in the 
planetary boundary layer. Both types of work serve not only 
to allow more realistic simulations but also prepare us to 
handle real Seasat data. 

The second extension is a set of experiments designed 
to determine the optimal assimilation period for Seasat wind 
data. Current practice is to assimilate data .at each model 
time step — every 10 m.inutes. Using longer assimilation periods 
is not only computationally more efficient but may yield better 
results because it reduces ^the extent of the boundaries of 
regions where data are" inserted. The model fields necessarily 
adjust incorrectly at a boundary between forecast and observed 
values. On the other "hand, too" long a period clearly degrades 
the data by increasing the difference between the time the 
observation was made and the time it is assimilated- Preliminary 
results thus far suggest that a 3-hour assimilation period 
yields the same results as a 10-minute one. Table 2 shows errors 
in forecasts obtained from initial conditions in which the data 
were assimilated in 1/2-hour, 3-hour, and 6-hour blocks. 
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Table 1 


Comparison of 3-day Forecast Surface Pressure Errors from 
Control (C ) , Error Free (SCM) and Random Error (SCME) Initial 
States (RMS in mb) . 
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TROPICAL SEA SURFACE TEMPERATURE: 

AN INTERACTIVE ONE -DIMENS lONAL ATMOSPHERE -OCEAN MODEL 

(E, Sarachik) 

Since the atmosphere and ocean are interacting systems, 
it is inappropriate to specify the sea surface temperature when 
dealing with the atmosphere or atmospheric anemometer level 
temperature and moisture when dealing with the ocean. All of 
these quantities should be determined interactively in terms of 
the external forcing: the solar ' constant . 

In the tropics, it is shown that the (cumulus) con- 
vective processes may be described by a one-dimensional cloud 
model. The near-surface ocean may similarly be described by a 
one-dimensional, mixed- layer model. The coupling is achieved 
through a sea surface flux budget combined .with the flux para- 
meterizations implied by the Monin-Obukhov similarity theory. 

This coupled one-dimensional atmosphere-ocean model 
is applied to the equilibrium situation wherein neither the 
ocean nor the atmosphere is heated or cooled. Since the ocean, 
lacking an internal heating or cooling mechanism, can only be 
heated or cooled through sensible heat fluxes through the sea 
surface, in equilibrium these fluxes must vanish. The atmosphere 
however maintains a stable lapse rate by balancing cumulonimbus 
heating against infrared radiative cooling. All water precipitated 
from cumulonimbus clouds must have evaporated from the sea surface. 
It is shown that this equilibrium system is closed and determinable 
solely in terms of the solar constant. 

For various values of the solar constant, the sea surface 
temperature, the flux of latent and sensible heat from the surface, 
the height of the tropopause, mixed layer, and trade inversion 
layer, and generally the entire vertical structure of the tropical 
atmosphere and ocean can be determined. The equilibrium sea surface 
temperature is shown to be relatively insensitive to changes in the 
solar constant, additional solar flux being compensated mainly by 
additional evaporation. 
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(not to scale) . (mixed layer height) ~600 m, Zipj (trade inversion) •^2 km. 

Zip (tropopause) ~16.5 km. The subsidence induced by deep cumulus -clouds 
acts to hold down the trade inversion. 


A TWO-LAYER DYNAMIC AND THERMODYNAMIC MODEL OF THE UPPER OCEAN 

(M. Cane and P. Schopf) 

A new model to be used to simulate the thermodynamics 
and circulation of the upper ocean is currently being coded and 
tested at GLAS/GMSF. The primary design goal is that the model 
be capable of state-of-the art simulations of those processes 
that determine the sea surface temperature. The model has been 
kept simple so that it is inexpensive to run simulations over 
several years. A coarse grid version has been structured to be 
compatible with the GLAS atmospheric GCM, This version is 
global and is intended to be used in coupled atmosphere-ocean 
experiments with the GCM. A finer-grid version can be used in 
a stand-alone mode either as a global model or to do experiments 
on a regional scale. 

The vertical structure of the model is depicted in 
Figure 1, The upper layer, the mixed layer, is treated as a slab: 
temperature and current velocities are taken as constant throughout. 
This is consistent with the bulk mixed-layer models now widespread 
(see Niiler and Kraus, 1977). The assumption that the layer is 
indeed well mixed is a good one for temperature (see Figure 2) but 
less so for currents. The depth of the mixed layer is- variable 
and is determined by the thermodynamics, turbulent mixing, and 
divergence of fluid due to large-scale motions. 

The layer below the mixed layer is taken to extend to the 
permanent thermocline with the deep ocean approximated as being 
at rest (cf., e.g., Charney, 1955). In this second active layer, 
the currents are again assumed to be uniform with depth, but the 
temperature is now taken to be linear with depth. This is the 
simplest structure that allows one to model the effects of static 
stability on mixed-layer entrainment. in 'both active layers the 
dynamics are governed by the primitive equations. The mixed-layer 
depth and temperature are also affected by local mixed-layer 
processes. These are modeled by a bulk formulation very similar 
to that of Kraus and Niiler (1977). Thus, in the model, as in 
nature, the surface temperature is influenced by local surface 
heating and wind stirring, horizontal advection, and upwelling, 
(Cooling due to upwelling is simulated by first having surface 
layer divergence thin the layer, thereby bringing the cold water 
below nearer to the surface where turbulent mixing processes 
entrain it into the mixed layer.) 

The model will be used to study the relationship between 
local thermodynamic and large-scale dynamical influences on sea 
surface temperature. 
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T = T 2 y = 0 

Figure 1. Vertical structure of the model. 
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T(‘C) 7 9 11 13 

S("M 32 33 34 35 



Figure 2. Typical variation of water properties with depth, 
November 1962. + = temperature, X = salinity, in 
parts per thousand by weight; 0 = density anomaly; 

" = 1000 density-1 where density is in g/cm^ and is 

measured at atmospheric pressure. From Grant, H.L. , 
A. Moilliet, and W. M. Vogel, 1968, Some observa- 
tions of the occurrence of turbulence in and above 
the thermocline. Journal of Fluid Mechanics , 34 , 
pp 443-448. 
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A DYNAMIC MODEL FOR THE OCEANIC MIXED LAYER: 

NON-ADVECTIVE VERSION 

(I. Fung and M. Cane) 

A necessary prelude to the coupled atmosphere-ocean 
general circulation model is an ocean model capable of reproducing, 
at least on the seasonal time scale, the sea surface temperature 
(SST) patterns given the atmospheric inputs of wind stress and 
net heat balance at the sea surface. This report outlines the 
first of a series of global models of the upper ocean currently 
under development and testing at GLAS/GMSF. 

The one-dimensional model of the mixed layer is fashioned 
after Niiler {1975). It basically consists of a mixed layer and 
a seasonal thermocline overlying a permanent thermocline. SST and 
mixed-layer depths are predicted. The mixed-layer equations in 
terms of the buoyancy content, total kinetic energy, potential 
energy and momentum of the layer are algebraic, and, hence, can be 
solved expeditiously. 

The one-dimensional model is applied to the GLAS 4° x 5° 
grid over the global oceans. Climatological wind stress (Hellerman, 
1967) and the net heat balance at the sea surface (Schutz and 
Gates, 1971) are used as forcing. The initial condition is given 
by the January NCAR climatological SST (Washington and Thiel, 

1970) and the January temperature profiles compiled by Miller 
(1974). The NCAR SST also serves as verification. The. time step 
used is 6 hours. It should be noted that integration for 1 year 
takes 5 minutes of CPU on the Amdahl 470. 

Figures 1 and 2 show the differences between the predicted 
and observed mean SSTs for January and March, respectively. They 
show that, although a 1-month integration (Figure 1) gives SSTs- 
that differ from climatological values by less than 1°C, a 3-month 
integration (Figure 2) gives systematic errors in the major current 
areas — where the effects of coastal and equatorial upwelling and 
advection by the western boundary currents are important. 
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A DYNAMIC MODEL FOR THE OCEANIC MIXED LAYER: ADVECTIVE VERSION 

(I. Fung and M. Cane) 

A 3-month integration of the one-dimensional model of 
the ocean mixed layer has shown systematic errors in the major 
current areas (Fung and Cane, 1978). This report outlines the 
second of a series of global upper ocean models that includes 
the effects of advection and uowelling. 

ih 

Denote the thermodynamic contribution to at, at 
by f) , g, respectively. Then the eguations governing the 
horizontally inhomogeneous 'mixed layer are 


3t 


V • (hu^) + h 


8hT 
m 

dt 


V • (hT u ) + q 
m~m ^ 


Since the currents are divergent, the magnitude and location of 
the upwelling velocity are calculated.. 

Because the mixing and advective processes are on 
different time scales, the splitting scheme (e.g., Marchuk, 1974) 
is appropriate. At each time step, the vertical mixing and 
horizontal current dynamics are thus treated separately. 

To use these equations, the currents have to be specified. 
A simple approach is to use a one-layer model where the mixed-layer 
velocity Uj^ is governed by the equation 

•It: ) = - V • (hu u ) + fkxhu ~ ~ f Vp dz + —(r - t ) 

3t ~m ~m~m -J ^ -~s ~e 

P -h P 


+ h /u - H(h) 6u) 


There is a hierarchy of approximations to the above equation. 
Only two will be listed here. 

1. The modified Ekman current given by 

1 

fkxhu_ = — (t - t ) 

~m - ~s ~e 
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where t = ku is provided at the bottom of the m-i-xed layer 
in 

so that the transport is <90® cum sole to the wind, as is observed. 
2. Modified Eknan and mean current 
Uj^ = u + 

v;here fkxhu^^^=~ t^=t - t 

- ~s ~e ~s ~s ~ 

P 


The mean current is calculated diagnostically from the density 
field. The level of no notion is taken at 1500 m; the current 
in the mixed layer is not very sensitive to its precise value. 

Experiments are underway to evaluate the inclusion of 
current dynamics on the sea surface temperature patterns'. 
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NUMERICAL SIMULATION OF THE GULF STREAM, GULF STREAM RINGS, 

AND MID-OCEAN EDDIES 

(T. N. Lee) 

The circulation and surface temperatures of the western 
North Atlantic were simulated with a primitive equation model to 
simulate the effects of mid-ocean scale eddies and their transports. 

Three important changes were made in the baroclinic 
model that was used by Semtner and Mintz (1977). First, a para- 
meterization of the subgrid scale diffusion of momentum and 
temperature was used, derived from the quasigeostrophic turbulence 
theory of Tokioka (1978). Second, a more realistic single-gyre 
wind stress forcing replaced the unrealistic 2 1/2-gyre pattern of 
wind stress that was used by Semtner and Mintz (1977). Third, 
because of the lack of observational data, friction at the lower 
boundary was neglected. 

The model has five levels and a horizontal grid size of 
37 km. The idealized model domain is a rectangular basin, 3000 km 
long, 2000 km wide and 4 km deep, oriented as shown in Figure la. 

The near shore side of the basin has a simple continental shelf 
and slope (shown in Figure lb), whereas the other sides are bounded 
by vertical walls. The model ocean is driven not only by a single- 
gyre pattern of steady zonal wind stress (shown in Figure la) but 
also by a Newtonian-type surface heating. 

In the experiment, the simulated Gulf Stream separates 
from the coast and continues eastward as an intense free jet. 

This Gulf Stream develops large amplitude transient meanders, 
some of which become cold-core cyclonic rings and warm-core anti- 
cyclonic rings that drift westward. Transient mesoscale eddies 
form in the broad westward-moving North Equatorial Current. The 
simulations of the flows and the surface temperatures are realistic 
in comparison with observations. 

Figure 2 shows the sea surface height at 40-day intervals. 
The contour interval is 20 cm. For the purpose of orientation, 
the geographical latitude lines at the northern and southern 
boundaries of the subtropical wind stress gyre are included, 
together with the 72 °W meridian-. These six sequence maps of the 
sea surface height demonstrate the life cycle of the Gulf Stream 
ring. On day 352, a meander of the Gulf Stream reaches sufficient 
amplitude to alter the flow from a generally eastward to a southerly 
direction, cold water is drawn south and a loop is formed. The 
loop moves westward and its area expands on day 392. On day 432, 
the ring is pinched off from the Gulf Stream and moves west of the 
72° meridian. On days 472, 512, and 552, the ring gradually decays, 
moves southwestward , and finally is entrained into the warm Gulf 
Stream. The ring has a life expectancy of approximately 300 days. 
The generation, maintenance, and decay of the rings in the model 
are comparable to observations. Newton (1961) has calculated that 
the heat exchanges across the mean Gulf Stream associated with the 
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Figure la. Boundary of the model domain superimposed on. a schematic 
picture of mesoscale variability in the western portion 
of the north Atlantic (Joint' Polymode Working Groups, 1976). 
Also shown in the zonal wind stress which drive the model 
ocean. 



distance FROM NORTHWEST COASTLINE (km) 

Figure lb. Model topography and grid box placement 
near the northwest coastline. 


118 



119 




120 






121 



Figure 2c. Surface heights fields (contour interval 20 cm) at 40-day intervals — day 432. 
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Figure 2d. Surface heights fields (contour 



20 cm) . at 40-day intervals — day 472. 
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Figure 2e. Surface heights fields ’ (contour interval 20 cm) at 40-day intervals — day 512. 
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separation of one ring from the main cold water mass, can be about 
200 to 400 X 10^® g-cal. Five to ten pairs of rings annually in 
the North Atlantic would be capable of transporting the required 
heat poleward by the North Atlantic Ocean (Oort and Vonder Haar, 
1976). Therefore, the correct simulation of the rings is very 
important to determine the heat transport by the ocean. 

This model appears capable of simulating the clima- 
tological properties of the Gulf Stream, its meanders, cut-off 
rings, and mid-ocean eddies. This is a step toward the eventual 
construction of a global oceanic general circulation model 
which, coupled to a global atmospheric general circulation model, 
can be used for interactive atmosphere-ocean climate variation 
simulations to help in understanding the global atmospheric-oceanic 
heat budget. 
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FORCED BAROCLINIC OCEAN MOTIONS III: 

THE LINEAR EQUATORIAL BASIN CASE 

(M, Cane and E. Sarachik) 

Previous work on the linear spin-up of an equatorial 
ocean is extended to include the specific effects of the north- 
south extent of the basin, thus allowing a detailed comparison 
of the analytic spin-up theory with numerical calculations. 

North-south modes in a &-plane equatorial basin are 
solved for both numerically and analytically. Simple approxi- 
mations are developed and the modes are classified by the location 
of the turning points relative to the zonal wall’s. The modes are 
illustrated for three cases: 

■ a. A symmetric basin whose zonal walls are 

distant from the equator (compared to the 
equatorial radius of deformation) ; 

b. A symmetric basin whose zonal walls are 
close to the equator; and 

c. An asymmetric basin, one wall of which is 
close to the equator and one far. 

Spin-up in response to x-independent winds in each of 
. these basins is then calculated in terms of four elements: the 

x-independent response to the wind; the eastern boundary response 
consisting of Rossby and Rossby-Kelvin modes needed to bring the 
u field to zero; the western boundary response consisting of 
Kelvin waves needed to bring the u field to zero; and the western 
boundary layer consisting of trapped short Rossby waves needed ■ 
to bring u to zero and meridionally redistribute mass. These 
elements are combined to describe how the steady (Sverdrup) 
solutions are approached. Special attention is paid to the fast 
planetary response: the equatorially confined Kelvin mode, the 

exponential anti-Kelvin mode with maximum amplitude on the zonal 
walls, and the Rossby-Kelvin mode, the latter two being modes 
not present on a meridionally open basin. 

Finally, a numerical model is run to illustrate spin-up 
for various winds in each of the three (a, b, and c) basins. 
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(d) \ u(1 ) 




The response to a uniform westerly wind of a zonally- 
unbounded ocean of meridional extent; 

(a) y = -5 to y = +5 

(b) y = -5 to y = +1.7 

(c) y = -1.7 to y = +1.7. 

The influence of the zonal costs is seen to be largely 
confined to a region within a local radius of 
deformation of the coast. 





Figure 2 . 


The first seven meridional eigenfunctions in a basin 
bounded by zonal walls at y = -5 and y = +1.7. 


has -been, non-c 


by the 


equatorial radius of deformation.) Note that the 
structure is severely modified by the near wall. 
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IV. CLIMATE 



SOME ASPECTS OP HYDROLOGY IN THE GLAS GCM 


(R. Godbole) 

Numerous simulation experiments have been performed 
using the GLAS GCM for winter ( January-February) and summer 
(July-August) months. The results of some of these experiments 
have been recently discussed by Charney (1978), Halem (1978), 

Herman (1978), and Shukla (1978). This paper is mainly concerned 
with the performance of the GLAS GCM with respect to handling 
some of the hydrological parameters. Specifically, the distri- 
butions of precipitation and evaporation are examined to see how 
well they agree with each other and with observations, and in what 
manner, if any, they are associated with other related parameters 
such as cloudiness and vertical velocity. 

For the purpose of the present discussion, the results 
of one summer month, namely July, were selected. Further, the 
paper's attention is restricted to the zonally-averaged fields only. 
This would mask the important regional features, but, nonetheless, 
it offers a convenient and powerful tool in understanding the 
persistent features of any quantity in relation to those of 
other related quantities. A detailed description dealing with 
the regional distributions of precipitation, evaporation, 
sea-level pressure,, etc., as obtained by the model for winter 
and summer seasons, and their comparison with observations has 
been documents'^ by Godbole (1978). 

Figure 1 shows the distribution of July precipitation 
(mm/day) as computed by the model. For comparison, the observed 
precipitation for the three summer months of June, July, and 
August compiled by Schutz and Gates (1972) is also reproduced 
in the figure. The agreement betv/een the two curves is fairly 
good. Both the profiles show a primary maximum due to the ITC2 
and two secondary maxima, one in each hemisphere, apparently 
associated with the frontal systems. The noteworthy feature is 
the dip in the precipitation amount in the winter hemisphere at 
a latitude of about 20°. The dip is more pronounced in the 
computed results. There is a secondary minimum in the observed 
as well as computed precipitation pattern in the summer hemisphere 
at about 40°. These precipitation minima are obviously due to the 
existence of the subtropical anticyclones which are .characterized 
by relatively clear skies arid large-scale subsidence. It can be 
inferred from the figure that the subsidence in the winter hemisphere 
is stronger than that in the summer hemisphere. This is known to 
be the case mainly because of the fact that the subsidence is 
associated with the descending branch of the Hadley cell, which is 
always strong in the winter hemisphere. This is clearly borne 
out from Figure 2 in which a zonally-averaged vertical velocity 
for the whole column is plotted along with precipitation and 
cloudiness as a function of latitude. All the quantities are 
model generated. By convection, negative vertical velocity 
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Figure 1. Zonally-averaged mean monthly precipitation 
(mm/day) computed by the model (thin line) 
and observed (thick line) . 



LATITUDE (DEG) 


Figure 2. Zonally-averaged mean monthly precipitation 
(mm/day) , cloudiness (percentage of time) 
and vertical velocity (10“^ mb/sec) as 
computed by the model. 
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represents upward motion. Cloudiness is expressed as a percentage 
of time any type of cloud is present in a vertical column. It 
is seen from the figure that between 55°N and 40°S all the three 
curves follow each other very closely, meaning that the upward 
motion is associated with relatively more cloudiness and more 
precipitation and vice versa. Over the polar latitudes, this 
conventional relationship apparently does not hold: specifically, 

the relationship between cloudiness and precipitation is reversed; 
more cloudiness being associated with less precipitation. Since 
the model physics precludes the existence of nonprecipitating 
clouds, this seeming paradox may be explained on the basis of 
saturation conditions and the actual availability of moisture. 

Over the polar latitudes, the almost saturation conditions would 
favor frequent formation of clouds, but the realization of 
moisture in the form of precipitation would be insignificantly 
snail. 

An important feature- to note from Figure 2 is that 
over the South Pole, the cloudiness is maximum . It immediately 
suggests that the South Pole should lose heat due to large 
outgoing longwave radiation, which is really not the case {Schutz 
and Gates, 1974). This may, therefore, well be the cause for 
unusually low temperatures generated by the model at the upper 
levels near the South Pole (Halem, 1978; Godbole, 1978). A 
furfher'analysis shown in Figure 3 supports this conjecture, 
because the two uppermost layers of the model (i.e., centered 
at 65 mb and 175 mb) exhibit a large amount of cloudiness near 
the South Pole. As a matter of fact, the longwave radiational 
cooling in the model (not shown) does show an increase in the 
cooling rate at the South Pole. 

Figure 4 shows the two major types of precipitation, 
namely, supersaturation precipitation and precipitation due to 
penetrative convection, which together contribute to the total 
precipitation. It is noticed that in the summer hemisphere the 
penetrative convection dominates the supersaturation condensation, 
whereas in the -winter hemisphere, the supersaturation processes 
overtake the convective processes. This is what should be 
expected. Supersaturation condensation may be considered to be, 
in part, due to frontal systems (besides due to ITCZ), which 
are more active in the winter hemisphere due to strong baro- 
clinicity. On the other hand, in the summer hemisphere, the 
atmosphere is convectively unstable ,' especially over the 
tropical belt where the precipitation mainly comes from the 
Asian and African monsoons. 

Figure 5 shows the computed surface evaporation. The 
corresponding observed evaporation derived by Schutz and Gates 
(1972) is also reproduced in the figure for the sake of comparison. 
It may be noted that the observed evaporation is not what is 
actually observed by any direct measurements; rather it is derived 
indirectly on the basis of the knowledge of other hydrological data. 
At any rate, it is found that the computed values of evaporation are 
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Figure 3. Zonally-averaged. distribution of cloudiness 
in a layer centered at 175 mb (dashed line) 
and at 65 mb (thick line) . 



LATITUDE (DEG) 


Figure 4. Distribution of supersaturation precipitation 
(dashed line) and penetrativ:e convective 
precipitation (thin line) in relation to the 
total precipitation (thick line) . The units 
are mm/day. 
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Figure 5. Zonally-averaged mean monthly evaporation 
(mm/day) computed by the model (thin line) 
and observed (thick line) . 


Dmparable with the observed values and the agreement between the 
wo curves is fairly good except that, in the Southern Hemisphere, 
he computed evaporation rate is somewhat underestimated. 

Finally, a brief comment is necessary on the hydrological 
alance of the model, which is represented in Table 1. The global 
anthly mean evaporation is found to be 2.69 mm/day, whereas the 
orresponding figure for precipitation is 2.34- mm/day. That means 
or a typical summer month of July, the global mean evaporation 
xceeds precipitation by about 13 percent. Th.is prompted GLAS to 
xamine further the daily behavior of global mean evaporation, 
recipitation, and specific humidity, which are the three major 
omponents of the hydrological balance. It was found that evapo- 
ation, precipitation, and specific humidity all decreased from 
he beginning until the end of the month with very little daily 
scillation in each of these. For convenience, the values on the 
irst and the last day of the month are shown in the table. It is 
een that evaporation exceeded precipitation all through the month, 
ut its rate of excess decreased during the month by about 4 percent 
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The corresponding decrease in the moisture content of the atmosphere 
is about 3 percent. These values,, however, must be viewed with 
some caution because of the Inherent computational errors and the 
•daily variation of a parameter itself. Moreover, certain parameters 
like ground wetness and precipitation runoff have not been taken 
into consideration. Nonetheless, the results of Table 1 do seem 
to point out that the model tends to maintain the hydrological 
balance within- a seasonal degree of accuracy. 


Table 1. Global Mean Evaporation, Precipitation 
and Specific Humidity for July Computed 
by the Model. 


Global 

Mean 

Evaporation 
(E) mm/day 

Precipitation 
(P) mm/day 

Excess 
(E)over (P) 
Percent 

Specific 

Humidity 

gm/kgm 

Monthly 

Mean 

2.69 

2.34 

13.0 


- 

‘ 

1st of 
July 

2.91 

2.49 

14.4 

4% 

decrease 

2.57 

3% 

decrease 

End of 
July 

2.49 

2.23 

10 .4 

of excess 
(E) over (P) 

2.50 

of specific 
humidity 
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RADIATION PARAMETRRIZATION AMD INFLUENCE OP SYSTEMATIC RADIATION 

DIFFERENCES ON THE GLAS GCM 

{M. L. Wu and L. D. Kaplan) 

Changes nade in the calculation of radiative heating can 
change the static stability of a model atmosphere and (instan- 
taneously) the temperature of the surface, v/hich in turn affects 
sensible heat flux and evaporation. Evaporation provides moisture 
into the atmosphere, which is the source of heat of condensation. 
Sensible heat flux alters the static stability of the atmosphere. 

The static stability affects the vertical motions field, which 
produces cloudiness and precipitation. The release of precipi- 
tation is accompan'ied by the release of latent heat, which in 
turn affects the motion field. Cloudiness again alters the 
radiative heating in the atmosphere and' surface temperature; 
this in turn affects the transport of latent and sensible heat 
into the atmosphere across the lower boundary. Because of these 
complicated interaction processes, systematic radiation differences 
v/ould be expected to be important in the modeling of the general 
circulation. A more realistic thermal radiation parameterization 
was, therefore, developed for the GLAS GCM, resulting in marked 
improvement in the model performance, as indicated by the results 
of impact studies summarized below. 

A thermal radiation routine was parameterized to replace 
an old radiation routine used in the GISS GCM. The parameterization 
is based on the statistical band model with a strong line version of 
the Curtis-Godson approximation. It includes the effect of water 
dimerization in the 8-13yn window region, and assumes a continuous 
distribution of specific humidity in the vertical. Calculations 
with this parameterization were found to be in good agreement with 
the theoretical line-by-line calculation of flux and flux divergence 
and to be within the observational error of long-wave radiative 
flux measurements. 

Winter and summer 10-week GCM runs were carried out with 
old and new longwave radiation parameterization, starting with real 
data for December 20, 1972 and June 18, 19,73, respectively, and 
the mean - .tistics for. the 30 days were compared for each of the 
two seasons. Figure 1 shows, for the initial condition, the 
longwave heating rates calculated with the old and new radiation 
parameterizations . 

The large difference in heating rates in the lower tropo- 
sphere at low latitude is due to inclusion of the waterdimer 
absorption in the 8-13;im region. The differences in the upper 
atmosphere are mainly due to different treatment of absorption in 
the 15ym carbon dioxide band. The new radiation parameterization 
also produces much less downward flux (approximately 30-50 percent) 
at the surface. 
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The introduction of the new radiation parameterization 
resulted in significant improvement in the model statistics and 
showed that the GLAS GCM is capable of simulating important clima- 
tological features remarkably well if a realistic radiation 
parameterization is used. 

This can perhaps be most easily seen by examination of 
the mass field. Figure 2 shows the mean sea level pressure 
distribution for the last 30 days of the winter run as obtained 
.with the old and new radiation, and includes for comparison the 
observed mean sea level pressure distribution as given by Taljaard, 
et al . , (1969) and. Crutcher and Meserve ( 1970). Figure 3 shows 
a similar comparison with the mean height field at. 200 mb level 
obtained by Crutcher and Davis (1969). 

In the winter experiments, the Icelandic low, the 
Aleutian low, and the mid 70 ceanic ridges in the Atlantic and 
Pacific are reasonably simulated in the new radiation run, 
especially the Icelandic low, v/hich is practically missing with 
the old radiation parameterization. The subtropical highs in 
the Southern Hemisphere are better simulated in the modified 
experiment. The westward tilt of troughs and ridges are also 
well simulated in the model with the new radiation. The sub- 
tropical highs in both hemispheres were simulated better in the 
■summer runs with the new radiation. 

The zonally-averaqed temperature field with the new 
radiation parameterization is closer to the observed temperature 
field in the tropics, but improvement is still required in the 
middle and high latitudes. At high latitudes, the temperature 
is too cold by as much as 20 °C at upper levels and by 10 °C at 
lower levels. Generally speaking, the temperature structure of 
the model atmosphere with the new radiation parameterization is 
more stable than that with the old radiation. The Hadley and 
Ferrel cells are stronger with the new modified radiation. 

A rainfall analysis, which was carried out cooperatively 
with Dr. Ramesh Godbole, showed the following differences between 
the old and new models. 

In the winter experiment, the modified model is able to 
simulate the precipitation associated with the ITCZ more realis- 
tically in regard to the horizontal extent and the amount than 
that of the old model. The wet areas in the northern Atlantic 
and Pacific Oceans coincide with major storm tracks in the modified 
model, but not so in the old model. On the whole, the precipi- 
tation rate in the new model over the oceanic regions is closer 
to what has been derived, from satellite observations (Rao, et al . , 
1976). The major dry areas off the west coasts of South America, 
Africa, and Australia are reproduced better in the modified 
model . 
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Figure 2a. 


Observed sea level pressure (mb-1000) for January after 
Jaljaard et al. (1969) and Crutcher and Meserve (1970) . 
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Figure 2b. Simulated sea level pressure (mb-1000), for February - Control 
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Figure 2c. Simulated sea level pressure (mb-1000) for February - 
Experiment , 
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Figure 3a. Observed 200 mb geopotential height (mb) for December 
January, and February after Crutcher and Davis (1969) 
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Figure 3b. Simulated 200 mb geopotential height (mb) for 
February - Control . 



Figure 3c. Simulated 200--rab geopotential height (mb) for 
February ’ - Experiment . 







In the summer rainfall, the significant feature is the 
precipitation associated with monsoon systems. The Asian SW. 
monsoonal rainfall is reproduced in the modified model much 
better than it is in the old' model. 

The eddy available potential energy (EAPE) increases by 
32 percent (10 percent) and zonal available potential energy 
(ZAPE) increases by 7 percent (27 percent) with the use of the 
new radiation parameterization in winter (summer) simulations. 

More detailed analysis of the energetics is in progress^, especially 
in the generation of EAPE and ZAPE. There are a few cases of 
worsening with the new radiation model, such that the Siberian 
high is worse with regard to intensity and location. On the whole, 
the worsening is not comparable with the large improvement 
mentioned in the previous discussions. 

These experiments have shown that accuracy in the 
radiation calculation is important in modeling the general 
circulation. Equally important is the computation time. A 
radiation parameterization that is accurate and fast in calcu- 
lation is highly desirable for the purpose of general circulation 
studies. GLAS is attempting to use a least squares approach to 
fit the emissivity to a polynomial. 

GLAS has also tested the role of ozone on the thermal and 
dynamical structure of the GLAS' GCM. The experiments were integrated 
for 42 days in a winter simulation. The last 30-day mean shows that 
the temperature in the upper layer (10-120 mb) in the tropical and 
subtropical regions decreases by about 5®-6°C when ozone-heating in 
the infrared absorption region is turned off and decreases by about 
8®-13°C when ozone heating in both the infrared and solar absorption 
region is turned off. This research is still in progress. 

The above experiments call radiation every 2 1/2 hours. 

If it is called every 1/2 hour, the zonally-averaged temperature 
field does not change much^ but a small though noticeable improvement 
is' obtained in the other fields, A detailed analysis of this 
interesting result will be carried out. 
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RADIATION-CLOUD FEEDBACK SENSITIVITY STUDY 


(Y. Sud and J. Shukla) 

A study was undertaken to examine the sensitivity of 
the radiation-cloud feedback mechanism using the general circu- 
lation model of the Goddard Modeling' and Simulation Facility.' 

The purpose of this study was to understand the importance and 
need for a time-dependent radiation-cloud interaction in a GCM. 

The GLAS GCM (previously the GISS GCM) has been described by 
Somerville, et al. (1974), Stone, et al. (1977), and more recently 
by Halem, et al . (1978). The latest version of the model, as 
described by Halem, was used in this study. An experiment for a 
90-day simulation was available from the in-house library of runs 
for the months of June, July, and August, starting from the 
initial data for May 15, 1974. 

The cloud parameterization in the standard model is 
interactive with the radiation. Clouds of supersaturation and 
convective types are generated in the model. The convective 
clouds are of three types: low levels 7 and 8, mid levels 5 

and 6, or deep convective between four specified levels 4 
through 9, while the supersaturation clouds occur if the rela- 
tive humidity exceeds the prespecified limiting value, usually 
100 percent. The short wave radiation is affected by clouds 
because the cloud reflects the incuident radiation, and multiple 
scattering within the cloud causes local warming due to the 
increased photon path. Long wave radiation, on the other hand, 
treats all types of clouds as black bodies for both absorption 
and transmittance. Typically, the presence of a cloud causes 
warming beneath the cloud, due to trapped long wave radiation. 

Features of the GLAS model used in this study that 
. differed from the earlier version described by Somerville, et al . 
(1974) and Stone, et al . (1977) were as follows: 

1. This model used a split grid. The standard 4°x5° 
grid changes to 4°xl0° and 4°x20°, respectively, at 66® and 82® 
poleward from the equator. The use of this grid spacing allowed, 
a 10-minut_‘ time step instead of a 5-minute time step. 

2. • The model used center differencing, MATSUNO-scheme 
instead of the TASU-MATSUNO scheme in the time-marching calculation 
in hydrodynamics. Besides being particularly suited for a split 
grid, this scheme eliminated an implicit horizontal diffusion in 
the TASU-MATSUNO scheme. However, a Shapiro filter on the tempera- 
ture and pressure fields was used to eliminate checkerboarding. 

3. The present model has fully interactive soil- 
hydrology. A simple though reasonable time-dependent ground 
moisture and ground snow and ice cover model was used. Snow 
and- ice albedos were calculated from the snow and ice cover 
depth on the ground. 
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4. The radiation in the earlier version of the model 
was replaced with the Wu-Kaplan radiation. This radiation is 
briefly described by Wu , et al . (1978). On the basis of identical 
summer and winter simulations with the old -and new radiations, 

Wu, et al . claimed the new radiation performs better than the 
previous radiation. 

The third month of the simulation period, i.e., August 
1974, was chosen for this experiment. The data set for fixed 
clouds was produced from the cloud output of the control run. 

First, a mean monthly cloudiness was obtained and then grid cells 
with highest percentage cloudiness were picked up to contain the 
clouds. Thus, clouds were allotted to a sufficient number of 
cells to conserve the zonal cloudiness separately over land and 
ocean and separately for each cloud type at each sigma level. 

The .new sensitivity run was then begun from August 1, 1974, using 
the data' set , boundary and initial conditions used in the control 
run. The only difference in the new run was an overriding state- 
ment that made radiation see the clouds from a fixed-cloud data set 
rather than those produced by the model for the purposes of pre- 
cipitation and convection. A mass-weighted cloudiness map for 
the control and experiment is shown in Figures la and- lb. The 
clouds produced from the dynamics of the experiment are shown 
in Figure Ic. Some of the significant results from this study 
are detailed in the following paragraphs. 

The sea-level pressure (SLP) fields for the mean monthly 
August looked quite similar (Figures 2a and 2b). In the difference 
fields (Figure 2c), it is noticed that in the Northern Hemisphere 
the new run gave SLP 4 mb deeper in the Pacific and 8 mb less 
intense in the Atlantic. In the Southern Hemisphere, there are 
three major cells with differences of -8, +4, and +8 mb, respec- 
tively, over Pacific, Atlantic, and Indian Oceans. 

The differences at the centers of these systems are 
larger than the model, variability (expressed as .standard deviation) 
obtained from a set of- four July predictability runs shown in 
Figure 2d. Nonetheless, there is no evident spatial correlation 
between locations of maximum cloudiness differences and maximum 
changes in SLP. 

m 

similar results were noticed in geopotential height 
fields. Notice differences ranging from -64 through +128 gpm in 
the Southern Hemisphere and -30 through +160 gpm in the Northern 
Hemisphere in Figure 3a. These differences are larger than the 
model variability from four July predictability runs shown in 
Figure 3b. However, positive height differences are consistent 
with positive SLP differences and vice versa. Again, there is no 
spatial correlation between locations of maximum changes in 
cloudiness and maximum changes in geopotential heights. 

The atmospheric energetics fields from these two runs 
were also examined. Zonally-averaged fields of evaporation. 
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Figure lb. Time- and mass -weighted fixed 
clouds in the experiment. 
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Figure Ic. Time- and mass-weighted clouds produced 
in the experiment. 
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(Experiment - Control) Aug 74 


SEA LEVEL PRESSURE (DIFF) 



Figure 2c. Differences in the sea level pressures 
(experiment - control) . 


SEA LEVEL PRESSURE ( standard Deviation) July Predictions 



Figure 2d. Model variability from four summer runs 
expressed as standard deviation among 
the four predictibility runs. 
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precipitation, temperature, specific humidity, and winds were 
compared. No major impact was found, therefore, these fields 
are not shown. 

In the two-dimensional fields, the most conspicuous 
differences are in the net outgoing fields of radiation. Obviously, 
it is caused by the cloudiness fixed at 100 percent at certain grid 
cells, which gave an inhomogeneous radiative transfer field. The 
fields, from control and experiment are shown in Figures 4a and 4b, 
respectively. Notice Figure 4a is a typical structure seen in 
observations (Halem et al . , 1978). 

The difference fields (experiment-control) of total 
evaporation and precipitation are shown in Figures 5a and 5b. 

Thick digits are for positive differences and thin ones are for 
negative differences. The only significant impact seems to be 
an increase in the evaporation in the western Pacific Ocean. 

Further investigation is necessary for this peculiar difference. 

In conclusion, it may be stated that the mean monthly 
sea level pressure and geopotential height fields do not show any 
systematic difference in association with changes in cloudiness. 

In certain areas the changes in these fields are larger than the 
natural variability of the model. As expected, the structure of 
the net radiation balance at the top of the atmosphere is very 
different in two runs. It is noisy in the fixed-clouds run. 

However, the significant changes are in the evaporation . fields 
in the western Pacific, These need further investigation to 
ascertain the interrelationship of cloudiness and evaporation. 
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RADIATION BALANCE (Control) Aug, 74 



Figure 4a. Radiation balance at the top of the atmosphere 
in (cal cm“2 min“^) for control. 


RADIATION BALANCE (Experiment) Aug. 74 



Figure 4b. Radiation balance at che top of the atmosphere 
in (cal cm“2 min“^) for experiment. 


152 



TOTAL EVAPORATION 


(Experiment - Control) Aug 74 



Figure 5a. Difference in the total evaporation in, 0 . l*itnn/day . 
Bold print is for positive difference; Tight 
print is for negative difference. 


TOTAL.PRECl PITATION (Experiment - Control) Aug. 74 



Figure 5b. Difference in the total precipitation in 0,l*min/day 
Bold print is for positive difference; light print 
is for negative difference. 








References 


Halem, M. , J. Shukla, et al , ^ 1978: Comparisons of observed 

seasonal climate features with a winter and summer numerical 
simulation produced with the GLAS general, c-i-rcul-ation- model. 
J.OC .Study Conf. on Climate Models, Washington, DC, April. 

Somerville, R. C. J. , P. Stone, M. Halem, et al . , 1974: The GISS 

model of the global atmosphere; J . Atmos. Sci . , 31 , 
pp. 84-117. 

Stone, P. , S. Chow, and J. Quirk, 1977: The July climate and a 

comparison of the January and July climate simulated by the 
GISS GCM. Mon . Wea . Rev . , 105 , pp. 170-194. 


Wu, M. L., L. D. Kaplan, and R. Godbole, 1978: Influence of 

systematic radiation differences on the dynamics of a model 
atmosphere. Third Conf. on Atmos. Rad., Davis, CA, Aroer. 
Meteor. Soc., June. 


154 



THE EFFECT OF CUMULUS FRICTION ON THE GLAS MODEL 
OF THE GENERAL CIRCULATION 

(M. He If and) 

A simple scheme to parameterize the vertical mixing 
of horizontal momentum has been added to the GLAS model of the 
general circulation. The scheme is similar to those of Ooyama 
(1971), Arakawa and Mintz (1974), Schneider and Lindzen (1976), 
and Anthes (1977). The cumulus friction term can be expressed as 

-^{Mc(U - yc)} , 

where Mc(p) is the total upward mass flux within all convective 
elements in the grid box at a given level, U(p) is the vertical 
profile of the ambient horizontal velocity, and Uc(p) is the 
in-cloud profile of velocity. The latter profile is determined 
by the assumption, 

^(McUc) = U , 

9p 3p ~ 

that horizontal momentum is conserved by a rising, entraining 
parcel' of air and by taking as a boundary condition 

Uc(P base) = U(P base). 

•Sj ^ 

Thus, in an entrainment layer the net effect of cumulus friction, 

^ {Mc(U - Uc)} = ' 

is that of the subsident advection of ^ which occurs between the 
convective elements. At the top of the cloud there is a detrainment 
layer, Ap, where the net cumulus drag is 


P top 

{Mc(U - Uc)\^P — {Mc(Uc - U)}p bottom 
9p ~ ~ ^ Ap ~ ~ 

p bottom 

When cumulus friction was added to the GLAS model, it 
had little effect in the tropical circulation. When it was added 
to the model in conjunction with an expanded cloud spectrum that 
included deep, ."hot tower" convective clouds, the effect was more 
dramatic. The January winter Hadley circulation accelerated by 
14 percent (compare Figures 1 and 2), the vertical profile of the 
mean meridional wind field became smoother (compare Figures 3 and 4) 
and the mean zonal shear of the tropics decreased (Figure 5). The 
Hadley cell strengthening due to cumulus friction is comparable in 


1 
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HEIGHT (mb) ‘ HEIGHT (mb) 



The January mean meridional stream function (x gm sec 

as simulated in the control. Negative regions are shaded. 
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Figure 2. The January zonal mean meridional wind field (m sec ) as 
simulated in the control. Negative (northerly) regions 
are shaded. 
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The January zonal mean profile of the 
zonal wind (m sec~^) as simulated in 
the experiment (solid line) and control 
(dashed line) and as observed by Oort and 
Rasmusson (1971) (dot-dashed line) . 






magnitude to tbe estimate given by Schneider and ^Lindzen (1977) 
and Schneider (1977) on the basis of an axisymmetric, steady- 
state model. Stone, et al . , (1974) also found a Hadley cell 
strengthening when they added a universal term for the vertical 
mixing of horizontal momentum to the GISS model, but they 
found that eddy kinetic energy in the middle latitudes was 
seriously damped. In the present experiment, cumulus friction 
actually caused a slight increase in eddy kinetic energy due 
to the development of a tighter meridional temperature gradient 
in the subtropical region. 

The strengthening of the Hadley circulation can be 
explained in terras of the atmosphere's angular momentum budget. 
Cumulus friction, through the subsidence and detrainment terms, 
carries relative angular momentum downshear from the poleward 
to the equatorward branch of the Hadley cell. The response to 
this imbalance is an enhanced Hadley gyre which now carries 
additional planetary angular momentum northward into the 
poleward branch and removes an additional amount of planetary 
angular momentum from the equatorward branch to the south. 

These results imply that the Hadley circulation is regulated 
by the atmosphere's angular momentum budget. The mean poleward 
temperature gradient, therefore, drives the Hadley circulation 
via the thermal wind relation. 
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REGIONAL COMPARISONS FOR WINTER AND SUMMER WITH THE GLAS GCM 

(M. Halera, J. Shukla, Y. Mintz, M. L. Wu, R. Godbole, 

G. Herman, and Y. Sud) 

Integrations were carried out for Northern Hemisphere 
winter and summer seasons using the GLAS general circulation 
model (GCM). Each integration was initialized from a real 
global analysis provided by the National Meteorological Center (NMC) . 

Initial conditions for the winter simulation were valid 
for OOZ January 1, 1975; the initial conditions for the summer 
simulation were valid for OOZ May 15, 1974. The data for the 
analyses were collected during the early Data Systems Tests (DST) 
conducted by NASA and made available to the Goddard Modeling and 
Simulation Facility by the Development Division of NMC. 

In the winter simulation, the Facility used the Posey 
and Clapp (1964) regional distributions of surface albedo and snow 
and ice cover. In the summer simulation, only three values of 
surface albedo (oceans [.07], vegetated land [.14], and deserts 
[.35]) were used in place of the Posey and Clapp data. 

The winter integration was carried forward for 60 days 
and the summer integration for 108 days. The simulated atmos- 
pheric state variables and diabatic processes were stored at 5-hour 
intervals to sample the evolving fields at every hour of the day 
(once every 5 days). The simulated seasonal climates described 
below are based on averages, sampled with the 5-hour frequency, 
over the last 30 days of each integration, namely February and 
August. 


Sea Level Pressure . Figures la and lb show the observed 
and simulated sea level pressure fields for February; Figures 2a 
and 2b show the respective fields for August. The observed fields 
are from Crutcher and Davis (1969) and are based on 100 years of 
observations. 

February. The Aleutian and Icelandic lows and the oceanic 
highs in the Northern and Southern Hemispheres are in good agreement 
with the observations. The position and depth of the Aleutian low 
is particularly good, as is the two-cell structure of the Icelandic 
low, although these are features which in nature show appreciable 
variations from year to year. The major shortcomings in the sea 
level pressure simulation are the excessively high pressure over 
the Arctic basin, the relatively weak and somewhat displaced Siberian 
high, and the failure to get the belt of closely spaced isobars 
between 40S and 60S. 

August. The features which are in good agreement are the 
high's over the subtropical oceans (which have hitherto been the 
most difficult feature to simulate in the Northern Hemisphere 
summer), the South Asia low and its extension across North Africa, 
and the low over the western U. S. The major shortcomings of the 
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Figure la. Observed sea level pressure (mb) for February, after 
Crutcher and Davis (1969) . 



figure lb. Simulated sea level pressure distribution (mb - 1000) 
for February. 
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simulation are that the Asian low is too deep and centered 
somewhat too far east, and the pressure is too high over 
Antarctica, 

Geopotential Heights . Figures 3a. and 3b show the 
observed and simulated geopotential heights of the 200 mb pressure 
surface for February; Figures 4a and 4b show the respective fields 
for August. The observed fields are taken from Crutcher and Davis 
(1969) . 


February. The features which are reasonably well simu- 
lated are the ridge .and trough located over the western and eastern 
coasts of North America, the ridge over the eastern Atlantic, and 
the strong Southern Hemisphere zonal flow. The major defect is the 
failure to get the tight packing of the height contours over Japan, 
which corresponds to the subtropical jet stream in that region, 

August. The features which are reasonably well simulated 
are the Tibetan high and its westward extension over North Africa, 
the ridge over Mexico, and the strong Southern Hemisphere zonal 
flow. The major ddficiency is that the observed ridge over the 
Gulf of Alaska is displaced to the center of North America. 

Surface winds . Figures 5a and 6a show the direction of 
the mean surface observed “wind for January and July; Figures 5b 
and 6b show the corresponding simulated fields for February and 
August. The observed fields are from Mintz and Dean (1952). 

These direct ahalyses of the surface wind field, as well as 
their upper- level counterparts, reveal the structure of the 
atmospheric circulation in the low latitudes much more clearly 
than can be inferred from the sea level and upper pressure fields. 

February. There is a good correspondence between the 
location of the intertropical convergence zone (ITCZ) in the 
simulated and observed surface wind fields. Particularly striking 
is the similar way in which the ITCZ varies its position with longi- 
tude: in each case having its main discrepancy in the ITCZ location 

is over the eastern part of the Pacific and over the Indian Ocean 
where the simulated ITCZ lies along the equator, whereas the observed 
ITCZ is a., about 7®N latitude over the eastern Pacific and between 
5° and 20°S latitude over the Indian Ocean, The flow of the surface 
air from the northeast and the" southeast toward the ITCZ (the "trade 
winds") is also well simulated. 

The centers and the asymptotes of anticyclonic outflow 
in the subtropical latitudes, which are the origins of the trade 
winds, are in good general agreement, espe.cially in the Southern 
(summer) Hemisphere, and the simulated and the observed outflow 
centers disa.gree to the extent that the simulated and observed 
surface pressure fields (Figures la, lb) disagree. Thus, the 
simulation shows a separate high pressure center and a separate 
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Figure 3b. Simulated 200 mb geopotential height (m) for February. 
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Figure 5a. 


Direction of the observed wind at the earth's surface for 
January, after Mints and Dean (1952) . 




Figure 5b. Direction of the simulated wind at the earth's surface 
for February. 
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Figure 6a. Direction' of the observed wind at the earth's surface for 
July, after Mintz and Dean (1952) . 
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Figure 6b. Direction of the simulated wind at the earth's surface 
for August. 
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center of anticyclonic outflow in the western north Pacific, 
which the long-period averaged observed pressure and wind 
fields do not show, but such a double-cell structure is 
observed in the north Pacific in some individual years. 

August. The change in the location of the ITCZ from 
February to August is very much alike in the observed and simu- 
lated surface wind fields. Particularly striking is the change 
over South America, central and east Africa, the Indian Ocean, 
south Asia, and the westerji Pacific. Again, it is the simulated 
summer hemisphere (now the"" Northern Hemisphere) in which the 
subtropical anticyclonic outflow is most clearly defined. Over 
the north Indian Ocean the surface winds have correctly changed 
from northeasterly to a southwesterly direction. 

Radiation Balance of the Earth-Atmosphere System . 

Global satellite observations of the geographical distribution of 
the radiation balance of the earth-atmosphere system are available 
for selected time intervals for a 10-year period. Such data can 
be used as diagnostic measures of the GCM ' s ability to simulate 
the prognostic variables that affect the radiation balance, namely, 
temperature, cloudiness, vapor distribution, and surface reflec- 
tivity. 

Table 1 describes the radiative properties of the clouds 
in the GLAS model. The short wave -radiation parameterization recog- 
nized eight types of model clouds to which it assigns the indicated 
optical properties shown in Table 1. The long wave radiation 
parameterization treats all of the clouds as "black." 

Figures 7a and 8a show the net radiation balance derived 
from Nimbus-3 satellite observations for the last 2 weeks in July 
1969; Figures 7b and 8b show the simulated radiation balances 
for February and August to be consistent with our other charts. 

(The solar declination angle changes by 8° in the second half of 
January and 10® in February.) 

February. The geographical distribution of the net 
monthly radiation balance that the model simulates reasonably 
well are the centers of maximum heating over the oceans in the 
low latitudes of the Southern Hemisphere and the relative minima 
over southern Africa. A ma>or difference, however, is the position 
of the line of zero heating, which in the simulation is about 5® 
to 10® too far north. 

August. The model simulates the northward shift of the 
lines of zero heating and the latitude- of maximum heating. 

Especially noteworthy is the fact that the model produces a minimum 
in the summertime radiational heating over north Africa. This 
causes the Laplacian of heating in the horizontal direction and 
its effect on the horizontal circulation, identical to the observed 
field. This holds even though the sign of the heating at the center 
of that region is not the same as in the observations. Also well 
simulated by the model is the minimum of the radiational heating 
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Table 1. Radiative Properties of Clouds in the GLAS Model. 


Cloud Type 

Layers 

Analogy 

Optical 

Thickness 

Albedo 

Super saturation 

2 

Ci 

1 

.12 

Super saturation 

3 

Ci 

2 

-21 

Supersaturation 

4 ■ 

As 

4 

.34 

Super saturation 

5 or 6 

As 

6 

.44 

Super saturation 

7 or 8 or. 9 

St 

8 

.51 

Penetrating Convection 

4-7 or 5-8 or 6-9 

Cb 

32 

. 81 

Mid-Level Convection 

5 or 6 

Ac 

8 

.51 

Low-Level Convection 

7 or 8 

Cu 

16 

.68 


NOTE: For long wave heating ' rates , all clouds are "black." 


( 1 ) 


Layers are '•110 mb thick, and layer one corresponds to interval between 10 mb and 120 mb. 








Figure 7a. 


Observed radiation balance (cal cm“2 min~^) of the 
earth-atmosphere system during the period January 21 to 
February 3, 1970, after Raschke, et al. (1973). 
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Figure 7b. 


Simulated radiation balance (cal cm min"'-^) of the 
earth-atmosphere system for February. 








over northwestern North America, The model does not reproduce 
the observed heating minimum off the west coast of North America, 
This may be due to the fact that the model does not simulate the 
planetary boundary layer stratus clouds. 

Global Radiation Budget , Table 2 compares the observed 
and simulated globally-averaged net radiation balance and its 
separate components for February and August and also shows the 
respective mean monthly planetary albedos. The observed fields 
are from Raschke, et al , ^ (1973),' It can be seen that the simu- 
lated albedo is about 10 percent higher than the observed. The 
incoming solar radiation differences are small and are probably 
mainly the result of using a solar constant of 1.95 cal cm“^min“^ 
for the model simulation. 

The simulated absorbed solar radiation shown in Table 2 
is in close agreement with observation, accounting for the good 
agreement in the planetary albedo. Table 2 shows that the simulated 
outgoing long wave radiation is about 12 percent less in winter and 
17 percent less in summer than observed. As a result, the net 
radiation balance is higher in the simulation by .037 cal cm~^min~^ 
in winter and .055 in summer. Several reasons may be responsible 
for this. First, clouds are assumed to be uniformly black to long 
wave radiation over the entire grid area (400 x 400 km and 100 mb) 
regardless of the type'" of clouds. Second, the model generates 
excessive high-level cloud cover in the equatorial regions, 
contributing not only to a higher albedo but also lowering the 
outgoing long wave radiation; the outgoing long wave radiation 
responds more strongly to an overestimate of high-level clouds 
than does the short wave radiation. 

Although the net radiation balance is positive in both 
seasons, the global mean temperatures of the simulation are con- 
stant to within 0.2 percent. From Table 3, it is seen that the 
net imbalance at the surface of the earth-atmosphere system is 
comparable with the imbalance of radiation at the top of the 
atmosphere as shown in Table 2. The storage in the atmosphere 
and over land is negligible in both seasons (not shown in Table 3). 
This implies that the oceans with their prescribed temperatures, 
and melting of permanent sea-ice, can act as a sink in the model 
to absorb the excess radiation balance calculated at the top of 
the atmosphere. 

Global Precipitation . Observations of the geographical 
distribution of rainfall over the oceans were recently acquired 
from the microwave instrument. Electrically Scanning Microwave 
Radiometer (ESMR) , on the Nimbus-5 satellite. A rainfall atlas 
derived from these data has been prepared by Rao, et al . (1976) 
with 4°x5° resolution. Figures 9a and 9b show the observed and 
model-simulated fields of rainfall for February; Figures 10a and 
10b show the respective fields for August. 

February. The general features of the monthly mean 
rainfall over the oceans are well simulated by the model. The 
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Table 2. Global Radiation Budget of the Earth- Atmosphere System (cal cm min ). 



Solar Radiation 

Albedo 

Long Wave Radiation 

Net 


Incoming 

Absorbed 


Outgoing 

Radiation Balance 

Observed 
(Jan. 21 - 
Feb. 3, 1970) 

.501 

.361 

.283 

.337 . 

.024 

Simulated 
(Feb . ) 

.513 

.356 

. 305 

.295 

.061 

Observed 
(Jul. 16 - 
31, 1969) 

.472 

.339 

.281 

.354 

-.015 

Simulated 
(Aug . ) 

CO 

00 

* 

.334 

. 315 

.294 

.040 
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Table 3. .Heat Balance at the Surface of the Earth (cal cm ^ min~^) 


February 

August 


Net Flux 
Incoming Solar 
Radiation 

Net Flux 
Upward Long 
Wave Radiation 

Net Flux 
Upward Sensible 
Heat Flux 

■ Net Atmospheric 
Heating by ' 
Condensation 

Net 

Surface 

Balance 

.278 

-078 

.040 

.101 

.059 

.251 

.073 

> .048 

.094 

» 

.036 


This is approximated by net evaporation 
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Figure 9a. ESMR-derived oceanic rainfall rate (.1 min/hr) for February 
1975^ after Rao, et al . (1976). 
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Figure 9b. Simulated oceanic rainfall rate (.1 ram/hr) for February. 
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Figure 10a. ESMR-derived oceanic rainfall rate (.1 mm/hr) for August 
1974, after Rao, et al . (1976) . 
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figure 10b. Simulated' oceanic rainfall rate (.1 irnn/hr) for August. 
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simulated field shows the bands of precipitation along the ITCZ 
in the Pacific and the Indian Ocean, although somewhat too far 
south in the. central Pacific. The model, however, is deficient 
in its simulation of the Atlantic ITCZ rainfall. Over the North 
Pacific, North Atlantic, tropi'cal Indian Ocean, and middle lati- 
tude Southern Hemisphere oceans, the simulation shows larger 
precipitation rates than the observations. 

August. The general features of the simulated fields that 
agree well with observations are the equatorial band of precipitation 
across the Pacific, the precipitation in the Bay of Bengal, and the 
intense precipitation over the tropical and subtropical western 
North Pacific, although the simulated amounts are too large. Again, 
the model failed to produce any significant precipitation along 
the Atlantic ITCZ. Another discrepancy is that the model produced 
a band of precipitation oriented NW-SE over the ocean east of 
Australia, which does not show in the observations. The model 
produced a band of rainfall over all of the middle latitude Southern 
Hemisphere oceans that was not observed. 

Conclusion . Monthly mean statistics of some atmospheric 
fields, as obtained from two seasonal numerical simulations with the 
current version of the GLAS GCM, were compared with the correspond- 
ing observed climate statistics. Although the simulated fields are 
for a single realization (with climatological mean boundary condi- 
tions) and the observations are for a number of years, the Facility, 
nevertheless, has attempted to infer from these comparisons the 
ability of this model to simulate seasonal climate. 

The geographical distributions of the simulated and 
observed fields were in good agreement for most of the large-scale 
features of sea level pressure, geopotential heights, horizontal 
velocities, radiation balances, and precipitation. The zonally- 
averaged statistics of winds and temperatures and the poleward 
transports of zonal momentum and sensible heat for the most part 
also agreed with the observations. It is recognized that there 
are also some significant quantitative differences between the 
simulated and observed fields. 

Moreover, many of these deficiencies are discussed in a 
more comprehensive paper. Comparisons of Observed Seasonal Climate 
Features with a Winter and ^nmer Numerical Simulation Produced with 
the GLAS General Circulation Model, which was included in the 
Proceedings of the JOC Study Conference on Climate Models (April 3-7, 
1978, Washington, D. C.). Future seasonal simulation studies will 
address the geographical distributions of the transient eddies or 
variances in time, space-time spectral analysis, and simulated 
model statistics over comparable ensembles of seasons with appro- 
priate boundary condition variations. 
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SENSITIVITY OF THE GENERAL CIRCULATION TO ARCTIC SEA ICE BOUNDARIES 

(G. F. Herman and W. t. Johnson) 

A series of numerical experiments was conducted with 
the Goddard general circulation model to test the. hypothesis that 
year-to-year variations of wintertime Arctic Sea ice boundaries 
affect hemispheric circulation features. Many observational 
studies into possible climatic effects of sea ice have been 
carried out since the beginning of this century, and these are 
summarized in detail by Herman and Johnson (1978). The predominant 
effect of ice cover variation is to alter the exchange of heat and 
moisture between the atmosphere and the surface of the ocean, which 
is due to the strong insolating property of sea ice. There is an 
additional effect on the radiation ' balance due to the larger albedo 
of ice as compared with ocean, but it is of secondary importance 
during the winter because of small values of insolation at high 
latitudes. 


Large excursions about the average wintertime sea ice 
cover occur in the East Greenland Sea, Barents Sea, Davis Straits, 
sea of Okhotsk, and Bering Sea. To test the effect of ice margin 
variations, GLAS performed a series of experiments in which the 
only control parameter was the total ice cover in these north 
Atlantic and north Pacific seas. The differences between the 
respective maximum- and minimum ice conditions then represent 
precisely the effect of ice boundary conditions on the GCM 
calculations. 

The maximum and minimum sea- ice extent in the north 
Atlantic is estimated based on ice charts from the British Climato- 
logical Service, which date from 1961 to 1977. Ice limits in the 
north Pacific were obtained for the 1973-1977 period from the U.S. 
Navy Fleet Weather Facility. These ice boundaries are shown in 
Figures 1 and 2, along with the Goddard GCM gridpoint representati^on . 
The area enclosed by the broken line represents an envelope of 
extreme ice conditions. (In general, ice limits do not vary 
synchronously in all Arctic seas.) 

Except for minor modifications, the version of the model 
that was used is essentially that described by Somerville, et al . 
(1974) and Stone, et al . (1977). Model initial conditions were 
obtained from the OOZ NMC analysis for January 1, 1975. A single 
realization of a January-February climate is simulated by inte- 
grating the model from January I initial conditions until 
February 14. Model output at 12-hour intervals (OZ and 12Z) is 
averaged for the period January 14-February 12. 

The ensemble mean control (i.e., minimum ice conditions) 
is estimated as the arithmetic average of six January-February 
simulations, and the mean anomaly (i.e., maximum ice conditions) 
is similarly estimated from two simulations. The individual 
simulations that comprise the control and anomalies differ from 
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one another by small, random initial values. These initial state 
differences are amplified by the barpclinic, barotropic, and 
convective instabilities of the model, and result in different 
mean monthly fields. (These instability mechanisms are some of 
the processes that cause one January to differ from another in 
the real atmosphere.) 

The difference in sea level pressure (minimum minus 
maximum ice) is shown in Figure 3. V^hen sea ice is at its minimum 
extent in the Barents Sea or Davis Straits, sea level pressure is 
lower than when at maximum extent. In general, the stronger 
surface heating over open ocean is associated with intensified 
cyclonic activity or, alternatively, more low level convergence 
and upper-level divergence, and lower sea level pressure. Higher 
pressure is found between Greenland and the British Isles when 
there is less ice present, and this change is of opposite sign 
to the Davis Straits and Barents changes. This results from mass-- 
conservation: surface divergence and high pressure southeast of 

Iceland accompanies the convergence and lower pressure to the 
east and to the. west of Iceland. Sea level pressure is lower in 
the Sea of Okhotsk with less ice, and a compensating pressure 
rise occurs in the Gulf of Alaska. 

The statistical significance of sea level pressure 
changes is estimated by comparing differences to the "inherent 
variability" or "noise" of the model. The high latitudes of the 
Northern Hemisphere, especially the regions of the north Atlantic 
and north Pacific where ice margin variability is the largest, 
have extremely high noise levels, i.e., standard deviation of 
sea level pressure are in excess of 5 mb. When normalized by the 
standard deviation of the sea level pressure (Figure 4), regions 
where the " signal-to-noise" ratio (difference divided by standard 
deviation) exceed two are found in the Gulf of Alaska, Sea of 
Okhotsk', Davis Straits, Barents Sea, and Icelandic north Atlantic. 
Large changes also occur in the subtropical north Atlantic and 
north Pacific. When averaged over latitude bands, sea level 
pressure changes are significant at the 6-percent level between 
50-70 °N, and at the 11-percent level between 10-30®N. 

urge changes also occur in the mid-troposphere . At 
500 mb the effect of maximum sea ice is to cause 2°C colder 
temperatures in the northern north Pacific, in the north Atlantic 
between Greenland and the British Isles, and over northern and 
central North America. These temperature changes correspond to 
decreased geopotential height surface at 500 mb -of over 100 gpm. 
These differences are three times greater than the standard 
deviation of the 500 mb temperature field in the -north Pacific, 
and twice as great over North America. The 500 mb geopotential 
differences are significant at the 3-percent level betv/een 50-70 ®N 

Changes of opposite tendency in the pressure and geo- 
potential fields in eastern subtropical Atlantic and Pacific 
accompanied the changes in the Icelandic and Aleutian regions. 
These changes may imply that changes in the mass field in high 
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Figure 1. Approximate maximum and minimum sea ice limits for North Atlantic 
seas, 1961-1977. Solid line is Goddard GCM grid representation 
of envelope of wintertime ice margin extremes. 



183 



Figure 2. Approximate maximum and minimum sea ice limits for North Pacific 
seas, 1973-1977. 
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latitudes and subtropical regions are interrelated, or "tele- 
connected" (cf., Namias, 1958). On the basis of these model 
results it is concluded that ice margin anomalies are capable of 
altering local climates in certain regions of the high- and m-id- 
la-titudes-/ and are also capable of exerting synoptic influence 
on a hemispheric scale. 
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SPECTRAL AND ‘SPATIAL ENERGETICS 


(J. Tenenbaum) 

GLAS has carried out calculations of the interannual 
variation of nodel and observational energetics and studied the 
effects of alternative calculations of (o, the vertical pressure 
velocity, on these results. Both of these efforts are directed 
toward the goals of: 

• Extracting the maxinum possible information for 
predictions in the range of 2 to 20 days 

o Specifying in greater detail where the model is 
failing to reproduce the observed behavior of 
the atmosphere. 

It is known that there is general agreement between 
the model and the atmosphere, particularly in the integrated ' 
values of mean and eddy potential and kinetic energies. (Most 
general circulation models have difficulty with the eddy kinetic 
energy.) The GLAS method uses spectral techniques in one and 
two dimensions to pinpoint the areas of disagreement. In 
particular GLAS is examining closely the long waves that would 
seem to be least affected by resolution and have the most 
kinetic energy. 

Studies of the interannual variation shed light on 
the model and observational variances, allowing detailed traces 
by wave number of where energy is being dissipated and generated. 
Three pairs of model and observation-runs have been analyzed 
over the Northern Hemisphere (see Figure 1). There is some 
suggestion that the model has difficulty generating eddy kinetic 
energy, Ke,' but can more readily follow dissipation of Kg. 

GLAS is currently analyzing runs aimed at verifying this conclu- 
sion by studying other years and using a model with a more 
realistic radiation code. 

The primary source of eddy kinetic energy is the 
conversion from eddy available potential energy? unfortunately, 
it is directly proportional to w and its spectral transitions, 
both of which are poorly determined by mars convergence in 
the observational motion field. 

Such an approach suffers from the errors inherent in 
the observational wind field. In many situations the signal- 
is entirely buried in the noise. The diagnostic W equation 
substitutes the solution- of an elliptic partial differential 
equation involving the mass field to eliminate some of these 
errors , 
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Preliminary results from this calculation are the following 
(see Table 1). The conversion from eddy potential to eddy 
kinetic energy is moderately sensitive to these errors at 
the lowest wave numbers, and quite so beyond wave number 4. 
As expected, the diagnostic omega equation significantly 
smooths the resulting . omega fields, perhaps, unrealistically 
so. GLAS is currently investigating the latter possibility. 
Inclusion of moisture has little effect, though without a- 
definite direction. 


Table 1. 


Conversion of Eddy Potential to Eddy Kinetic Energy 
(units w/m^).. Data' for OZ, December 20, 1972. 







Wave Number 





Case 

■ 

2 ' 

3 

4 

5 

6 

7 

8 
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10 

Mass 

Convergence 

0.74 

0.53 

0.65 

■ 0.37 

0.19 

0.24 

0.15 

0.18 

-0.03 

0.43 

w Equation, 

- No- Moisture 

1.23' 

0,04 

-0.25 

0.15 

0.02 . 

-a.OL 

O.Ql 

t0.03 

-0.06 

0.06 

to Equation, 
With Moisture 

1.24 

0.48 

0.51 

0.48 

0.02 

0.00 

0.04 

0.08 

0.01 

0.06 
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EFFECT OF A PACIFIC SEA SURFACE TEMPERATURE ANOMALY ON THE 
CIRCULATION OVER NORTH AMERICA 

(B. Bangaru and J. Shukla) 

During the fall and winter of 1976-77, sea surface 
temperature (SST) in the north Pacific was characterized by 
abnormally cold temperatures in the central and western portions 
of the north Pacific with a warm pool located off the west coast - 
of the U. S. (Figure 1). Namias (1978) has suggested that the nor-th 
Pacific SST anomalies may have been one of the multiple causes 
of the abnormally cold temperatures in eastern North America 
during the 1976-77 winter. This study attempted to test this 
hypothesis by conducting a numerical experiment with the GLAS 
general circulation model. 

The model used in the present study is the one described 
earlier by Somerville et al . (1974) and Stone et al . .(1977). 

A modified version of the model has also been presented in the 
proceedings by Halem et al . (1978). The experiments were per- 
formed by first integrating the GLAS GCM for 45 days with .the 
climatological mean SST and the observed initial conditions valid 
for January 1, 1975. This integration is referred to as the 
Control run (C). The climatological SST field was then changed 
by adding a time invariant anomaly field and therefore, although 
the climatological SST varied with season, the imposed anomaly 
field remained constant with time. The structure of the imposed 
SST anomaly was the same as shown in Figure 1, but the numerical 
values were exaggerated by considering them in “C rather than °F. 

The model was integrated again for 45 days; this run will be 
referred to as the Anomaly run (A). 

Control Run C and Anomaly Run A were further repeated 
with their respective boundary conditions in SST, but the initial 
conditions in u, v, T, and surface pressure were randomly perturbed. 
The spatial variation of the random perturbations corresponded to 
Gaussian distributions with zero means and standard deviations of 
1°C in temperature, 4 m/s in horizontal wind components and 3 mb 
in surface pressure over land points and 2°C in temperature, 8 m/s 
in horizontal wind components, and 6 mb in surface pressure over 
ocean points. 

These runs are referred to as Initial Condition Perturba- 
tion Runs Cj and A^ for Control (C) and Anomaly (A), respectively. 
Four additional Initial Condition Perturbation Runs were already 
made by Spar et al . (1978) which are referred to as C 2 r C 3 , C^, and 
C^. The atmospheric parameters for the six runs are referred to as 
C, Cj, C 2 r C^, C 4 , C 5 , as the mean control (Cm) and the average of 
two anomaly runs A and A^^ as mean anomaly (Am). In order to assess 
the response of the imposed SST anomaly on the atmospheric circula- 
tion, the differences between the mean control and the mean anomaly 
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Figure 1. Observed sea surface temperature anomaly 
(°F) during January 1977. 
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runs averaged for the period between day 15 through day 45 are 
examined and the magnitude of the differences is compared with 
respect to the standard deviation among the mean monthly values 
for the six control runs for this period. 

Figure 2a shows the differences betv/een the mean anomaly 
and the mean control runs for the 700 mb temperature. The largest 
differences are found along the International Date Line, which is 
the longitude of maximum negative anomaly in SST, and the eastern 
parts of Canada, and the northern u. S. A positive anomaly is found 
over Greenland. Occurrence of colder 700 mb temperatures over the 
cold anomaly and warmer 700 mb temperatures over the warm anomaly 
can be explained by considering the .radiative , sensible and latent 
heat fluxes. A detailed examination of the model-generated heat 
fluxes showed that the evaporation and convective cloudiness 
decreased over the cold, anomaly. Reductions in the moisture and 
the sensible heat flux, which reduces the convective clouds, causes 
a reduction in the latent heat of condensation and this can also 
cause cooling of the atmospheric temperatures. Consistent with 
the decrease in the convective cloudiness, it was also found that 
the solar flux at the surface increased but the long wave flux at 
the surface decreased. Since the sea surface temperatures were 
prescribed, this feedback does not operate in the model. It 
should be noted, however, that the net effect of these diabatic 
heat sources and sinks is such that the heat source (sink), 
which is initially confined to the lower boundary, finally 
extends into the troposphere. 

The most interesting feature in this figure is the 
coldest 700 mb temperatures centered around SON, 75W. This is 
clear evidence of the dov7nstream response of the model to the SST 
anomalies in the Pacific. Figure 2b shoves the ratio of the 
differences and the standard deviations among the six control 
runs. These ratios are the signal- to- noise ratio,. Two areas of 
maximum signal to noise are found, one over the anomaly itself 
and the other centered around SON, 75W. The positive temperature 
anomalies over Greenland were found to have a signal- to-noise 
ratio less than two. 

In the earlier numerical experiments with the NCAR model 
(Kutzbach et al . [1977], Chervin et al . [1976]) no significant 
downstream effect was noticed. Although the two general circula- 
tion models are not identical (especially with regard to their 
simulation of the transient eddies), in GLAS ' opinion, the primary 
reason a downstream effect was observed in -the GLAS model was due 
to the nature of the spatial structure of the SST anomaly. 

The results of these numerical experiments support the 
hypothesis by Namias that the SST anomalies over the Pacific 
during the 1976-77 winter may be one of the contributing factors 
towards very cold temperatures over the northeast U. S. and warm 
temperatures over the Alaskan region. Since this was not a 
coupled ocean-atmosphere model, SST anomalies were assumed to 
persist for the whole period of integration. This assumption is 
justified by the observed fact that the SST anomalies do persist 
for a time period of several months. 
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These results show clear evidence of a downstream 
response of the model atmosphere to the imposed SST anomalies 
over the Pacific. It is interesting to note that for some 
levels the downstream effect away from the anomaly is more than 
the effect over the anomaly itse.Lf. Since the natural variability 
of the mid-latitude atmosphere is rather large, the standard 
deviation among six predictability runs was calculated and 
the ratios of the mean differences to the standard deviations 
indicated that the downstream effects are indeed -significant. 

The results of the present study suggest that the 
downstream orographic barriers may contribute to amplify the SST- 
generated atmospheric perturbations. In other words, an oro- 
graphically ' forced atmospheric flow regime may be more sensitive 
to further modification by the perturbations generated by the 
SST anomalies. This qualitative conjecture is suggested by those 
numerical experiments in which it is found that the effects of the ^ 
north Pacific SST anomalies are most prominent on the western and 
the eastern sides of the North . American continent. This is of 
special significance for the winter of 1976-77 because, as pointed 
out by Namias ( 1978), the anomalous flov; pattern over North America 
v/as in phase with the normal winter flow pattern. Therefore, this 
was interpreted as an inphase amplification of the normal winter 
circulation, which 'is orographically and thermally forced by the 
mountains and the mean heat sources and sinks ._ The purpose of 
this numerical experiment was not to simulate the actual events 
of 'the 1976-77 winter but only to examine the effects of similar 
SST anomalies on the model atmosphere. It should be noted that 
the atmospheric conditions in the mean control run corresponded 
to the winter of 1975. 

Detailed examination of the model-simulated daily ^charts 
indicated that in general the anomaly runs did not show a persistent 
blocking pattern near the west coast of -the U. S. as observed during 
the 1976-77 winter. However, there were periods of 3-5 day duration 
in which the ridge over the west coast of the U. S. was highly 
intensified and such events were more frequent in the anomaly run 
compared to the control run. Blocking situations do occur but they 
do not last too long. It is not clear whether this is a manifesta- 
tion of some model deficiency or it is because the .model integrations 
are not performed long enough for several interannual cycles. 
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LOW-FREQUENCY PLANETARY WAVES: SPACE-TIME SPECTRAL -ANALYSIS 

OF A WINTERTIME SIMULATION OF THE GLAS MODEL 

{D. M. Straus, J. Shukla, and W. Byerly) 

- Space-tine spectral analysis provides a useful diag- 
nostic tool, applicable both to observations and nodel output. 

It allows one to separate various space and tine scales, and 
thus distinguish motions believed to have different origins 
or distinct dynamics. One can study medium-scale (wave number 
5 to 10 ) "baroclinic" motions at mid-latitudes separately 
from the largest scale (planetary) motions (Pratt, 1975, 

Pratt and Wallace, 19 76-) and observe the dominant time scales 
involved with each. Further, it is possible to look at the 
change in the character of motions of particular tine and 
space scales with height (Deland, 1973), with latitude (Hayashi 
and Colder, 1977), or both. 

This paper presents preliminary result's of a space- 
time analysis of a 60-day (wintertime) simulation of the GLAS 
model. Various fields (geopotential height, geostrophic wind 
components) were spectrally analyzed with respect to zonal 
wave number and frequency at a number of latitudes and model 
levels. The focus here is on the low-frequency, planetary 
scale motions which, in mid-latitudes, dominate the total 
variance of, for instance, the geopotential height (Pratt, 

1975, 1977). 

Space-time spectral analysis, as discussed by Pratt 
(1976), differs from straight-'forward Fourier analysis, in 
longitude and time in that the method recognizes the stochastic 
character of the atmosphere. One expresses the total space- 
time variance of a field q(X^t) as 

[q 1 = 1 / dt(q(X,t))'^ = ^ ^ [E^g(k,+w)4Egg(k,-m] (1) 

4irf k=0 o)=0 

where [ ] indicates a zonal average, an overbar indicates a 
time average, and the total length of record being considered 
has been normalized to 2tt. EqCT(k,w) represents the contri- 
bution to the total variance by motions of wave number k and of 
a band of frequencies centered at to. Here (0>o refers to 
westward movement. The stochastic nature of the atmosphere 
necessitates a certain amount of averaging over frequency in the 
calculation of Eqg(k,o)j in order that the’ resulting spectrum 
be satistically stable. With conventional techniques, this leads 
to a large amount of smoothing (e.g., Pratt, 1975, who smooths 
over a frequency range of .125 cycles per day) . A technique that 
allows much finer frequency resolution is the maximum entropy 
method (MEM, Hayashi, 1977a, Ulrych and Bishop, 1975, Smylie, 
et al . , 1973). This method was used in this study. 
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Recent work (Pratt, 1976, Hayashi, 1977b) has 
enphasized the importance of distinguishing between the con- 
tributions of standing waves and those of propagating waves 
to E^_(k,o)). Large-scale standing waves, (presumably forced 
by mountains and continent-ocean heating differences) , have 
fixed nodes, and consist of coherent eastward" and v/estward 
propagating waves of equal amplitude and frequency. They have 
the property that (k,+u ) = £L (k,-o) ) . However, the equality 

of (k,±w) does not necessarily imply standing waves: 

Physidally distinct, incoherent westward and eastward moving 
v?aves of the same frequency and amplitude also yield ( k , -Ho ) - 

Eqg(k,-o)). Distingushing betv/een these situations retires the 
calculation of the coherence between eastv^ard and westward propa- 
gating waves (Hayashi, 1977b). Since this is very difficult to 
calculate with MEM (Hayashi, 1977a), it has not been carried out 
in this study. 

Follov/ing the above discussion, one .defines (dropping 
the subscript 'qq' ) : 

Total variance 


T(k,u) = E(k,(o) + E(k,-to) = E^(k,to) + Eg(k,co) , • to>0, (2) 

* 0 
which represents the contribution to [q ] from motions of 

wave number k and frequency w and -to. 

Propagating variance 

PR(k,w) = E(k,(o) - E(k,-io) = E^(k,io) - Eg(k,w) , w>o, (3) 


which represents the net tendency towards propagation’ of 
motions of v/ave number -k and frequency ±(jj{PR>0 refers to 
westward propagation) . 

The renaip.ing variance, T(k ,o)) - I PR( k , lo) I , consists in 
general of standing waves, incoherent traveling wave, and 
random fluctuations. Hayashi and Colder (1977) indicate, in 
particular, that standing waves do make an important contribution 
on the time and space scales of interest here. 

Results are presented for E(k,±o 3 ) for the geopotential 
height (({)) for the longest planetary waves (wave numbers 1-5) 
and frequencies corresponding to 60 day-10 day periods, at 
mid-latitudes and at several levels. 

Table 1 presents the tota'l space- time variance T(k,w) 
for various periods and wave numbers at 54N, 510 mb. The most 
striking feature is the band of high variance extending from - 
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Table 1. Total Variance T(k,(u) (Arbitrary Units) 54 n, 510 nb 


VJave No. 


Pe.r.iod (.Days) 



60 

30 

20 

15 

10 

1 

523,551 

199,485 

33,136 

20,742 

18,206 

2 

95,251 

.132,363 

8 6,10 7 

52,889 

■ 21,548 

3 

62,741 

6,174 

109,157 

199,838 

14,433 

4 

8,033 

18,527 

40,134 

80,200 

5,663 

5 

6,592 

5,723 

7,733 

6,659 

8,230 


Table 2. 

Propagating 

Variance PR(k 

, (jo) (Arbitrary Units) 

54N, 510 mb 

Wave No . 


Period 

(Days) 



60 

30 

20 

15 

10 

1 

+473,886 

+136,443 

- 3,168 

+ 2,714 

+ 4,571 

2 

+ 17,722 

- 52,423 

-53,526 

- 27,198 

+ 6,691 

3 

- 34,337 

- 27,227 

-67,366 

-191,333 

-12,866 

4 

- 2,940 

- 9,743 

-24,345 

- 75,836 

- 5,021 

5 

+ 1,192 

- 2,268 

- 6,164 

- 5,392 

- 6,450 
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low frequency, low wave number to higher frequency and wave number. 
This feature, also observed in the atmosphere (Pratt, 1975, 1977), 
indicates that the largest waves have distinctly preferred time 
scales associated v?ith them, tine scales of interest for long- 
range weather nrediction. 

The propagating variance PR is presented in Table 2 for 
various periods and v/ave numbers at 510 mb, 54N. Disturbances 
of wave number 1 in the 60- and 30-day category propagate pre- 
dominantly westward. This is true at most but not all of the mid- 
latitudes and model levels. This is observed in the atmosphere 
(Pratt and Wallace, 1976, Fraedrich and Bottger, 1978), but was 
not found by Hayashi and Colder. In the observations, the dominance 
of i)?estward propagation extends over wave numbers 1 to 3 and over 
all periods greater than 10 days. This general trend is not 
uniformly reproduced in Table 2, and this is so for other latitudes 
and pressure levels. Apparently,. eastv7ard~moving disturbances 
are somewhat over-represented in the model. 

It is useful to compare the vertical structure of the model 
disturbances to that of the atmosphere. One way of accomplishing 
this is to find evidence in the model data of the nodes deduced by 
Pratt and Wallace (PW) for the atmosphere. They find that much of 
the low-frequency (periods of 10-20 days)', low-wave number variance 
of mid-latitude geopotential and temperature can be explained in 
terms of two modes: In the first westward-moving mode, | (j) | is 

non-zero at the ground and increases with height. The second node 
is eastward moving and is characterized by having | (})| vanish at 
the ground and increase more sharply with height. Pratt and V7allace 
describe the first modes as external, equivalent barotropic, and 
the second mode as more baroclinic in nature. If these modes are 
important in the model, one expects the low levels to be dominated 
by westward motions, with the strength of the eastward motions 
increasing with height. 

In Table 3, E^(k,w) and Eg(k,w) are shown for several 
latitudes and heights. The values of Eg were averaged over 
wave numbers 1-3 and periods of 20, 15, and 10 days. It is clear 
that if Pratt's modes are present, they are partially obscured by 
other types of motion: E^^(k,ii)) does increase with height, but 

Eg(k,w) is still large at 730 mb. In aqreement with PW, the 
relative importance of, eastward-propagating waves does increase 
with height: 

Eg(175 mb)-Eg(730 mb)>E^(175 mb)-E^(730 mb) 

at all latitudes, and mb) > at 30“ and 54N. 

Ep(730 mb) E (730 mb) 

^ w 

The latter inequality does not hold at 70N due to the relatively 
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Table 3. Average Values of Eg (Arbitrary Units). 






Eg 


Mb 

3 ON 

54N '70N 

3 ON 

54M 

70N 

17 5 

7,443 

29,702 30,935 

16 ,364 

121,748 

61,647 

510 

2,849 

11,921 ■ 10,459 

8,596 

49,863 

10,084 

730 

2,575 

9,450 5,272 

4,239 

25,042 

11,573 

1 

Table 4. 

r 

Average Values of PR(k,u) 

(Arbitrary 

Units) . 



Mb - 

30N 

54N • 


70N 

• 

175 

-8,923 

-92,546 


-30,655 


510 

■ -5,441 

-37,942 


+ 3,872 


730 

-1,664 ■ 

■ -27,469 


- 6,301 
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large value of Eg at 7 ON, at 7 30 nb. In comparing Table 3 to 
the propagating wave nodes of PW, one is assuming that standing 
waves do not contribute, for they would have to be subtracted 
from Eg, E^, before the latter give information on propagating 
waves . 


An alternate assumption is that the standing waves 
contribute as much as' they can v/hile being consistent with the 
data. This means taking the standing wave strength to be Eg(E^) 
if E >Eg (Eg>E ) . Subtracting the "standing wave amplitude*^ 
2Eg(2E^) from rhe total variance T gives in this case the 
propagating variance PR. 

Values of PR(k,dj) for various latitudes and pressure 
levels, averaged as before., are presented in Table 4. Again, 
the variance of geopotential corresponds to predominantly 
eastward movement, even at low levels, and this tendency 
increases with height {except between 730 mb and 510 mb at‘70N). 

Interpreting the model data with respect to the tv;o 
modes of PW, one can conclude that the baroclinic mode (second 
mode) is too strongly represented compared to the equivalent 
barotropic mode (first mode). To more carefully interpret the 
model data, it is necessary to calculate the coherence between 
eastward and westward moving components, so that the role of 
standing waves can be understood. This calculation is presently 
underway. 


Space-time spectral analysis of a wintertime simulations 
of the GLAS model shows that the low-frequency, low-wave number 
disturbances behave realistically in some ways. Their total 
variance is relatively strong -in a band extending from low wave 
numbers and frequency to high wave numbers and frequency. In 
addition, there is evidence of westward propagation in the longest 
wave length, lowest frequency category, although this tendency is 
not as widespread as in the atmosphere. 

Comparing the vertical structure of the model to the 
two modes of PW believed to be important in the atmosphere, one 
finds that the model's baroclinic mode is too strong compared to 
its equivalent barotropic mode. 
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A COMPARISON OF GROUND HYDROLOGY PARA^^ETERI ZATIONS 
USING INPUT DATA FURNISHED BY A GCM SIMULATION 

(J. D. Liri/ J. J- Alfano, and P. Bock) 

Three ground hydrology parameterizations proposed for 
use in an atmospheric general circulation model (GCM) are 
described in this paper. Controlled experiments were conducted 
for an off-line ground hydrology model for these parameterizations 
using input data furnished by a 2-day GCM simulation. A com- 
parison of results indicates that the parameterizations have 
significant effects on the change of ground wetness and actual 
evapotranspiration for a wet, an average, and a dry cell located 
in the continental United States. 

In the existing atmospheric general circulation models, 
the ground hydrology has been developed for the purpose of 
evaluating the source terms in the governing equations which 
are generated from the moisture, heat, and momentum fluxes across 
the atmosphere-land interface. The hydrology that was first 
incorporated into a GCM by Manabe (1969) was based upon the 
parameterization: (1) actual evapotranspiration as a function 

of soil moisture and potential evaporation ,- and (2) a certain 
moisture storage capacity in a ground surface layer or plant 
root zone. The dynamics of moisture-heat transport processes 
and the variable characteristics of the soil-water-plant system 
are grossly oversimplified in this scheme. Nevertheless, this 
simple approach has been shown to be capable of producing the 
basic features of hydrology and climate on a global scale 
(Manabe and Holloway, 1975) and in semi-arid regions (Charney, 
et al . , 1977) and remains in use' today in the primary GCMs 
(Manabe, 1969; SomerVille, et al . , 1974; Washington and 
Williamson, 1977; Gates and Schlesinger, 1977). This study 
was an attempt to evaluate the sensitivity of parameterization 
related to soil moisture in a ground surface layer and actual 
evapotranspiration by off-line numerical experiments using 
input data furnished by a GCM simulation. The parameterizations 
chosen for the study are three variants of Manabe 's scheme, 
which have been used in the UCLA model (Arakawa, 1972; Arakawa 
and Mintz, 1974); the GISS model (Charney. et al . , 1977; Lin, 
et al . , 1978); and the UC parameterization, a revised version 
proposed for use in the GLAS GCM (Lin, et al . , 1978). 

In a rectangular cell of a ground surface layer, the 
key parameter used in the ground hydrology (Manabe, 1969) is 
the soil moisture in a ground surface layer which is available 
to satisfy the evaporative demand of the atmosphere. In the 
UCLA model, ground wetness (w') is defined as the soil moisture 
(w) in a fraction of the maximum available soil water (w^^) .in a 
ground surface layer of unspecified thickness. In terms of 
soil properties (Salter and Williams, 1959), then 
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w 


w 


e 


(la) 


- 0 


wp 


w 


m 


fc 


- 0 


WP 


and = p^Sh (0jj, ■* 9^) 

where p is the mass density of water; S the soil porosity; 
h the thickness of ground surface layer; and 0, 0^^, and 0 
are volumetric moisture content and its value at field capacity 
and permanent wilting point,^ respectively. In the GISS and UC 
parameterizations , ground wetness is taken as volumetric moisture 
_ content (0)' in a fraction of saturation so that 

w' = 0 ■ (,1b) 


and Wj^ = p^Sh. 


The variable soil properties can be characterized by varying the 
maximum available soil moisture in the UCLA model and the field 
capacity and permanent wilting point in -the GISS and UC 
parameterizations , as shown in the following table;- 


Table 1. Soil Properties. 



UCLA ' 


GISS 

— 


UC 


Soil 

•Type 

Wm (g) 

Wjn (g-) 

6fc 

0wp 

Wn, (g) 

®fc 

CD 

*P 

Loam 

10 

66.4 

.25 

.10 

50 

.25 

.10 

Sand 

10 

4.0 

.09 

.025 

25 

.09, 

.025 


With the definition of ground wetness in (la) and (lb), a moisture 
conservation equation may be derived for a ground surface layer as 

8 S 

^ " ®s ■ ^o “ f^w (2) 

9t ot 


where P is precipitation, E evapotranspiration, R_ runoff 
including surface runoff and percolation, and S^^ the storage 
of snow and ice. 
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The actual evapotranspirat ion in (2) is scaled from 
the potential evaporation by means of a scaling factor g such 
that 


^ Ep 


(3) 


and 


^p ” Pa'-D^’^a' ^*3* (Tg) q^] 


(4) 


in which p^, V^, and are the air mass density, horizontal 
velocity, and specific humidity of the atmosphere at a reference 
elevation; q*(T ) is the saturation specific humidity at the 
ground temperature T ; is the momentum transfer coefficient 
assumed, to be equal to the moisture transfer coefficient. In 
the UCLA model. 


B = 


1, 1/2 < w' < 1 

2w', 0 < w’ < 1/2 


(5a) 


while in the UC parameterization. 


3 = 


1 

e -.0 


, < 0 < 1 


wp 


0.P -0 

fc wp 


, 0 < 0 < 0 ^ 

/ wp — fc 


(5b) 


The version of parameterization used in the GISS model and being 
studied here is the following: 


6= (Bpiatit+ Ssoll* ^ (1- "°CK) (5c) 

where 

Bplant = 1- - exp {[-0.1 x VPM x exp (16.75 x VPM)/Ep] x PLANT} 
Bgoii = 1* “ exp [-0.01 X VSM X exp (16.75 x VSM)/Epl 
VPM = 0.1 X VSM/USM + 0.2 X [(VSM - USM)/ASM] 

VSM = 2. X 0 X 
USM = 6^ 

A-SM = 0fc - V 

ROCK = fraction of impervious surface (being ignored in this 
study) 
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I.t is seen that the scaling factor is a function of soil 
properties, potential evaporation and solar radiation. The 
PLANT factor depends on the solar angle Z as 

PLANT = min [1.0, cos{ Z) /O .2588] 


The runoff in (2) consists of the surface runoff and 
the horizontal and vertical percolation out of the cell-layer 
control volume, representing the "loss" of soil moisture and thus 
unavailable for evapotranspiration. The runoff for the UCLA model 
is given by 


= [p3 + (10(1 - - 10{1 -w’) (6a) 

for the UC parameterization 

Rq = P(0.1 + 0.96^*^) (6b) 

and for the GISS model 

Rq = P - 0.5 X BBFACT , P 2 . x BBFACT 

, . (6C) 

= 0.5 X P^/BBFACT , P < 2. X BBFACT 

in which the coefficient BBFACT is expressed in terms of the 
following parameters: 

BBFAC = 2. X FILTF X exp [-1.75 X (ADINP - 0.5)] 

FILTF = [1.77 X exp{-0.036 x (Tg - 273.2)} + 0.44562]'^ 
and ADINF = 1.25 x (2. x 0 -0.2) - 0.5 


The term accounts for the storage of snow and ice 

which will 'accumulate when the temperature drops below 273. 2°K 
and melt when the temperature exceeds 273.2°K. The amount of 
melt is calculated by the heat balance to keep the surface 
temperature at 273.2 °K. 

The ground temperature is intrinsically coupled with 
the ground wetness in the GGM model. The ground temperature, 
the average temperature of the -ground surface layer, is governed 
by 
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LE + S. 
S 1 


(7) 


C 


3T _ 

g - 

dt 


R. 


N 


F„ - 


in which C is the bulk heat capacity, R», the net radiation flux 
at the ground surface, F the sensible -heat flux, LE the latent 
heat flux, and the" latent heat transfer due to melting of ice 
or snow or due to freezing of water. A detailed treatment of 
(7) has been documented by Tsang and Karn (1973). 

An off-line ground hydrology model (GHM) , forced by the 
atmosphere without feedback, was developed for these parameteriza- 
tions. The GCM furnishing the atmospheric input data contains a 
known ground hydrology which may or may not be the same as those ' 
in the off-line GHM for controlled experiments. 

Controlled experiments of the off-line GHM were con- 
ducted for the three parameterizations described in the previous 
section using input data furnished by a 2-day GLAS GCM simula- 
tion for a wet cell (Upper New England) , an average cell (Southern 
Great Plains) and a dry cell (Southwest region) in the continental 
United States. The soil for the dry cell is characterized by the 
sand type while for the average and wet cell by the loam type. 

The ground hydrology -used in the GCM is a version of the UCLA 
model, which assumes a constant, seasonally-averaged ground 
wetness. The input data which are extacted from the GCM global 
output consist of the surface pressure, temperature, specific 
humidity, wind velocity and net radiation flux. Only the initial 
values of ground wetness and ground temperature were used in the 
calculation and the rest of ground wetness and ground temperature 
together with the actual evapotranspiration were considered as 
diagnostic variables. 

The results of the calculated diagnostic variables, 
ground wetness, ground temperature and actual evapotranspiration, 
in each cell are shown for the parameterizations in Figures 1, 2 
and 3. In the 2-day run, there was no precipitation in the dry 
and average cell, but a total amount of 1.5 cm of rain fell on 
the wet cell. The daily average value of ground wetness, potential 
evaporation, actual evapotranspiration and the scaling factor 
for evaporation are listed in Table 2. 

The soil moisture in Figure 1 is described by the 
ground wetness defined in (la) and (lb). There is a significant 
change of ground wetness in response to evapotranspiration and 
light rain in the UCLA model even for the 2-day run, however, 
the UC and GISS ground wetness, remains at a constant level. The 
high sensitivity of ground wetness change in the UCLA model might 
be attributed to the value of as given in Table 1 and the 
value of 8 as indicated by the average values given in Table 2. 

It is clear that the length of a 2-day experiment is not suffi- 
cient to evaluate the change of ground wetness in the UC and 
GISS parameterization as an effect of soil type. 
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Volumetl'ic Moisture Content in Fractions of Saturation 




Aug. 18 '75 [hr] Aug. 19, '75 


Figure 1. Comparison of ground wetness. 


Ground Wetness in Fractions of Maximum Available Water 









Table 2. Comparison of Daily Average Ground Wetness, Potential 
Evaporation, and Evapotran’spiration . 


i 

j 

[ 

f 

1 1 

[ 

r 

1 

w' 


{ cm/h ) 

I ' 

1 

Eg ( cpi/h ) 

3 = 


Cell 

[ I 

Model i 

' ] 

i D^y 

. 1 Day 2 

Day 

1 Day 2 

Day 1 

Day 2 

Day 1 

Day 2 

( 

UCLA 

. 10 

.07 

. 19 

.19 

. 035 

.025 

. 18 

.13 

Dry 

UC 

.0 4' 

.04 

. 17 

.16 

.043 

.039 

.25 

.24 


GISS 

.22 

. 22 

.31 

.24 

.006 

.005 

. 019 

.021 


UCLA ' 

.39 

.37 

.062 

,061 

.048 

.046 

.77 

.75 

Average 

UC 

. . 17 

. 16 

.082 

.068 

.042 

*.03 4 

.51 

.5 


GISS 

.26 

.26 

.19 

. 16 

. 012 

.010 

.06 

.06 


i 

UCLA ■ 

.65 

.75 

.02 

. 055 

.02 

.055 

1 

1 

Wet 

UC 

.21 

. 21 

. 025 

.0525 

. . 0165 

.034 

.66 

. 65 


GISS 

.32 

.32 

.015 

.08 

.0075 

.024 

.5 

.3 
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The ground temperature in Figure 2 for all cases exhibit 
the diurnal heating pattern. The UCLA and UC results show good 
agreement. The cause of higher ground temperature, as much as 
10 °C, for the dry and average cell in the GISS model is believed 
due to the underestimated values of actual evapotranspiration. 

The values of actual evapotranspiration in the UCLA and 
UC parameterization are in good agreement despite the fact that 
the UCLA mode.l tends to provide higher values of g for soil 
moisture between field capacity and permanent wilting point for 
a wet cell and unreasonably low values for a dry cell. In the 
GISS model, the actual evapotranspiration was underestimated 
while the potential evaporation was overestimated. This result 
implies that the nonlinear behavior of potential evaporation in 
the parameterization of g in (6c) is too strong. 

This comparative study on the basis of 2-day numerical 
experiments indicates the following: 

a. The change of ground wetness in the UCLA model is 
more sensitive under both evaporation and precipi- 
tation conditions than normally expected during 
■such a short duration. The ground wetness in 

the UC and GISS parameterization appears to be 
relatively stable, however, longer numerical 
experiments are needed to test its sensitivity. 

b. In the GISS model, there exists a strong coupling 
between the ground temperature and actual evapotrans 
piration as„a result of the strong non-linear 
dependence of actual evapotranspiration on potential 
evaporation in the parameterization. 

c. The parameterization based on a constant maximum 
available water in the UCLA model causes the actual 
evapotranspiration to be underestimated in a dry 
cell and overestimated in an average cell. The UC 
parameterization appears to be better adapted for 
soil type specifications while the ground wetness is 
more relevant to the definition of soil moisture 
used in soil sciences. 
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NUMERICAL METHODS FOR METEOROLOGY AND CLIMATOLOGY 
{E. Isaacson, D. Marchesin, and G. Zwas) 

It has been estimated that about 30 percent of the error 
in meteorological forecasts is attributable to the inaccuracy of ' 
currently used numerical schemes. In order to minimize this 
error, GLAS has developed and tested two new numerical methods. 
These schemes are particularly designed to perform long terra 
integrations for climatological purposes, since they conserve 
the global quantities of total mass, total energy, and total 
potential enstrophy of the flow. One of the schemes is implicit, 
while the other one is explicit. The greater inherent efficiency 
of the. implicit scheme for fine spatial grid calculations arises 
from the fact that its time step At may be an integer multiple 
of AX/U, and is not required to be less than AX/(C+U), as is 
the case for explicit schemes, where AX is the smallest spatial 
interval size, U is the maximum speed of propagation of the flow, 
while C is the much larger maximum speed of sound waves in the' 
flow. In other words, a much larger time step is permitted for 
implicit schemes, but in order to derive this benefit it is 
essential that the. amount of computer time used to solve the 
implicit equations at any time step be a small multiple of the 
time needed to make one time step with an explicit scheme. The 
application of these schemes indicates that the implicit scheme 
does require less computer time than does the explicit scherae 
for fine grid sizes. 

In the following section both the implicit and explicit 
methods are used to solve the shallow water equations for a single- 
layered atmosphere over a smooth earth with comparisons of the 
solutions and computer time. In the next section, wi.th orography 
present, a comparison is made of the solutions produced by the 
explicit scheme with and without the conservation of the global 
quantities. Here it is observed that the solution produced, while 
conserving the global quantities, is smoother and more accurate. 

In the last section, a brief description is given of the numerical 
schemes. 


- Application of Schemes. The technique for devising an 
implicit scheme and the method for modifying both implicit and 
explicit schemes to conserve the global quantities were applied 
to the solution of the nonlinear shallow water equations for a 
single-layered, incompressible fluid over a spherical earth. As 
initial conditions, three subtropical highs were centered sym- 
metrically in each hemisphere on latitudes ±30“, and at longitudes 
0°, 120 °E, 240 °E. The initial wind speed was found geostrophically 
from the initial height field. The initial hydrostatic "pressure" 
(i.e., density-height) rose from 1036 (corresponding to 8500 
meters) to 1046 at the .center of the highs. The maximum of the 
initially geostrophic velocity was' about 13 meters/sec. Contour 
levels of the initial "pressure" field were plotted on a stereographic 
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chart of the Northern Hemisphere along with arrows inc^icating the 
initial wind speed and direction (see Figure 1). Figure 2 shows, 
after 24 hours, the appropriately scaled height fields produced on 
one-third of the Northern Hemisphere with a fine spherical grid 
using 360° in both longitude and latitude. Figure 2a was produced 
128 

by the explicit scheme; Figure 2b was produced by the implicit 
scheme. Note how closely the "high" cells resemble each other 
after 24 hours. The left— most column of numbers represent 
longitudes, i.e., row 73 corresponds to longitude 360 ° . 

The right-most column corresponds to the latitude 128 

nearest the equator and the second column from the left 
corresponds to the latitudes nearest the North Pole. 

For the explicit scheme a time step of 3 minutes was 
used, while in the implicit scheme an interval of 15 minutes 
was used. Table 1 lists the computing times (CPU) required on 
a 60-bit machine, the CDC-6600, for each method with both the 
fine grid and a coarse spherical grid of 360/64° for latitude 
and longitude. 


Table .1. Comparison of Computing Times on CDC-6600 (in minutes). 




Type of Scheme 



Explicit 

Implicit 

Revised 
Impl icit 

Coarse Grid 

3 

9 

.4.5 

Fine Grid 

24 

36 

18 


Note that the time step of 15 minutes was used' in the 
implicit scheme for both the fine and coarse grids; whereas in 
the explicit scheme, 6 minutes was used for the coarse grid and 
3 minutes was used for the fine grid. Both schemes produced 
solutions with comparable accuracy. Furthermore, it is possible 
to revise the implicit scheme to take . advantage of the special form 
of the hydrodynamic equations to cut the computing time of the 
implicit scheme in half. This feature of the equations of motion 
is also present in the full hydrodynamic equations used for general 
circulation models. With this improvement the implicit scheme would 
irun in much less time than the explicit scheme for fine grids. 

Orography Effects. A major source of computational 
errors is found to be the way in which current numerical methods 
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Figure 1. Contour levels of the initial pressure field. 
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Figure 2a. Height field after 24 hours - explicit scheme. 


73 21 I 2 J 2 444566077 / 7 /// 77//66666677 

74 2312 J 2444556677777 / 7777766 b 6 b 6 fti . 

75 23123 '2 34455567 /- 7 / 77 « d 77 7 66555666 

76 1312323 3 - 4556667777 MrtB 87 766545555 

77 1312221334556677778 8 . 888776433444 

76 131222 .3 3344556677788088775421 122 

79 1 3 1 22 I 222 134556 7 778-8868 7 752 -1 -3 -3 -2 -2 

80 131121211223 -* 5667 / 88 a 887750 -3 -4 -4 -4 -- 

81 131121 1 1010123566770888764 -2 -5 -5 -5 -5 -5 

82 1 3 1 1 Z 0 I a -1 -1 -2 -2 -I 0 3 5 6 7 7 0 0 8 8 7 6 4 -3 -5 -6 -6 -6 -6 

83 1 3 1 1 2 0 0 -1 -2 -3 -3 -4 -3 -3 -1 2 5 5 7 7 8 8 8 7 6 3 -4 -6 -6 -6 -6 -o 

84 2 3 1 1 1 0 -1 -2 -3 -3 -4 -5 -5 -5 -4 -2 2 5 6 7 7 7 7 7 5 1 -5 -6 -6 -7 -7 -7 

05 2 3 1 1 1 -1 -1 -2 -4 -4 -5 -5 -6 -6 -5 -5 -3 3 5 6 7 7 7 6 5 -2 -5 -6 -7 -7 -7 -7 

86 23111 - 1 - 2 - 3-4 -5 -5 -6 -6 -6 -6 -6 -5 -2 3 5 6 7 6 6 3 -3 -6 -6 -7 -7 -7 -7 

07 23111 - 1 - 2 - 3-4 -5 -6 -6 -6 -7 -7 -6 -6 -5 -2 3 5 6 5 4 1 -4 -6 -7 -7 -7 -7 -7 

08 2 3 1 1 1 -1 -2 -3 -5 -5 -6 -6 -7 -7 - 7 ' -7 -7 -6 -5 -3 2 4 4 2 -3 -5 -6 -7 -7 -7 -7 -7 

89 22111 - 1 - 3 - 4-5 -5 -6 -7 -7 -7 -7 -7 -7 -7 -6 -5 -3 -I -I -2 -4 -6 -6 -7 -7 -7 -7 -7 

90 22101 - 1 - 3 - 4-5 -6 -6 -7 -7 -7 -7 -7 -7 -7 -6 -6 -5 -4 -4 -4 -5 -6 -6 -7 -7 -7 -7 -7 

91 22101 - 1 - 3 - 4-5 -6 -6 -7 -7 -7 -7 -7 -7 -7 -7 -6 -6 -5 -5 -5 -5 -6 -6 -7 -7 -7 -7 -7 

92 22101-1 - 3 ‘-4 -5 -6 -6 -7 -7 -7 -7 -7 -7 -7 -7 -6 -6 -6 -5 -S -5 -6 -6 -7 -7 -7 -7 -7 

93 22101 - 1 - 3 - 4-5 -6 -6 -7 -7 -7 -7 -7 -7 -7 -7 -6 -5 -6 -5 -5 -5 -6 -6 -5 -6 -7 -7 -7 

94 22101 - 1 - 3 - 4-5 -5 -6 -6 -7 -7 -7 -7 -7 -7 -6 -6 -6 -5 -5 -5 -5 -5 -6 -6 -6 -6 -7 -7 
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Figure 2b. Height field after 24 hours'- implicit scheme. 
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treat flow over and around high and extensive mountain ranges, 
e.q.. Rocky Mountains, Andes Mountains, Himalayan Mountains. In 
order to eliminate this difficulty, Arakawa has advocated using 
difference schemes that conserve the total potential enstrophy, 
total energy, and total mass. GLAS has found it possible to modify 
any difference scheme to make the new schem.e conse,rve these global 
quantities. In particular, the effect that this conservative 
modification has is shown when orography is introduced into the 
shallow water model previously described. That is, three identical 
mountain ridges were introduced, centered respectively along the 
three longitudes 60'’E, 180?E, and 300®E. The maximum height above 
sea level of the mountain ridge was 8500 meters at the equator. 

16 

The mountain height decreased to zero at the poles. The initial 
height above sea level of the fluid layer was the same as shown' 
in Figure 1. The initial velocity was again found from the 
geostrophic balance on the sphere. The fluid tended to flow 
around the ridges by veering to the pole. In Figure 3a, it is 
seen that the solution produced affer 24 hours without the 
conservation of the global quantities has larger velocities 
and steeper height gradients near the pole than does the 
solution shown in Figure 3b, which is produced with the con- 
servation of total mass, total energy, and total potential 
enstrophy. In Table 2, the relative errors of these global 
quantities are indicated, e.g./ 

total mass (t) - total mass (0) , 
total mass (0) 

after t=24 hours and after t=48 hours for -the explicit scheme. 

In the first column the larger relative deviations are, found 
without the use of the conservative modification method; whereas 
the smaller relative deviations in the second column are produced 
with the use of the conservative modification method. 

The explicit scheme is a so-called leapfrog scheme in 
which the solution is advanced from (t-At) and t to the time 
(t+At). The first order time derivatives are replaced by 
differences centered at time t, while the first order spatial 
derivatives in latitude and longitude are replaced by fourth 
order, five point, centered .difference expressions. At the 
latitudes closest to the pole, the values of the solution 
determined by this difference scheme are smoothed by using the 
fast Fourier transform. This program has been adapted from the 
ideas of Kreiss and Oliger (1973) and Williamson and Browning (1973) 
and is described in Stoker and Isaacson (1975). 

The implicit scheme is, roughly speaking, a Crank- 
Nicolson type -scheme in which the solution is advanced from time 
t to time ( t+At) . Here the fourth order accuracy is obtained 
by using a Bade rational fraction in the spatial difference 
operators. After clearing fractions, by multiplying all terms 
by the denominator, a compact three point spatial difference 
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Figure. 3a. Contour levels of pressure-, after 24 hours, found 
without the conservation of the total mass, total 
energy, and total potential enstrophy. 



Figure 3b. Contour levels of pressure after 24 hours, found 
with the conservation of the total mass , total 
energy , and total potential enstrophy. 
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expression is obtained. The resulting simultaneous difference 
equations are simplified by writing the spatial operator as the 
product of a longitudinal factor and a latitudinal factor. These 
factors have a block tri-diagonal matrix representation that is 
easily invertible. This idea was proposed by A. Harten based 
on the methods used by Beam and Warming (1976, 1977) in aero- 
dynamics. Considerable effort was needed to find an .efficient 
and stable factorization method. Finally, it was necessary to 
use a fourth-order Shapiro filter at each time step to maintain 
stability. 


Table 2. Relative Errors in Total Mass, Total Energy, and 
Total Enstrophy - Produced With and Without the 
Use of the Conservative Modification Method. 





Conservative 

Modification 

Quantity 



Without ^ 

With 


At 

24 

Hours 


Total Mass 



140*10"^ 

0.2*10”® 

Total Energy 



203-10“® 

0.5-10”® 

Total Potential 

Enstrophy 

5766-10”® 

10.0-10”® 


At 

48' 

Hours 


Total Mass 



235-10“® 

0.5-10”® 

Total Energy 



260-10”® 

1.1-10“® 

Total Potential 

Energy 

17 

',750-10”® 

16.0-10“® 


The conservative modification method is described in 
general terms in Isaacson (1977). Here the work involved in modifying 
the solution found at time (t+At) is of the order of the amount 
of calculation used in one time step of -the ex plic-it scheme . 
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